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CHAPTER  1 


INTRODUCTION 


The  objective  of  this  thesis  is  to  present  the  results  of  an  investigation  into  the 
problem  of  producing  a  discrete  stochastic  time  series  that  has  any  desired  bispectral 
characteristic.  In  essence  this  amounts  to  an  inversion  of  the  bispectrum;  however, 
the  resulting  time  series  is  not  guaranteed  to  be  unique  in  any  way  except  that  it  has 
the  specified  bispectrum. 

Typically  most  papers  dealing  with  the  bispectrum  or  its  applications  are 
interested  in  how  one  estimates  the  bispectrum  from  a  sampled  subset  of  the 
underlying  random  process  [e.g.,  Brillinger  and  Rosenblatt,  1967b;  Rao  and  Gabr, 
1984;  Raghuveer  and  Nikias,  1986;  Nikias  and  Raghuveer,  1987  ].  This  thesis  is 
interested  in  exactly  the  opposite  problem,  i.e.,  how  does  one  obtain  or  estimate  the 
times  series  knowing  only  its  bispectrum.  While  this  is  not  a  well  studied  problem, 
fortunately  much  can  be  done  by  a  straightfoward  generalization  from  the  power 
spectrum,  which  is  very  well  known.  In  addition  there  is  one  advantage  to  inverting 
the  bispectrum,  that  is  not  present  in  the  bispectral  estimation  problem.  It  is  that  one 
can  assume  to  have  the  true  bispectrum  and  not  just  an  estimate  of  it,  unlike  the 
reverse  problem  where  one  cannot  generally  assume  to  have  all  the  representations 
of  a  random  process. 
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The  inversion  can  be  accoiiq)lished  fc^owing  a  procedure  first  suggested  by 
Wolinsky  [1988].  The  starting  point  is  ro  assume  a  model  for  the  time  series.  There 
are  many  models  one  could  choose  for  this  purpose.  However,  most  involve  veiy 
complex  and  often  unsolvable  relations  between  the  model  parameters  and  the 
bispectrum.  But  since  the  choice  is  open  it  only  makes  sense  to  choose  one  that  will 
make  the  inversion  as  easy  as  possible.  Acctmlingly  the  model 

x(t)  -  ti(t)  +  e(t)  +  II  g(r,s)e(t+r)T)(t+s)  (1.1) 

has  been  chosen,  for  reasons  which  will  be  explained  more  fully  in  chapter  3.  Here 
x(t)  is  the  desired  time  series,  e(t)  and  nCO  are  Gaussian  independent  identically 
distributed  (i.i.d.)  random  variables  with  zero  mean  and  variance  o^,  and  g(r,s)  is 
called  the  kernel  function.  This  is  a  slighdy  different  model  than  the  one  suggested 
by  Wolinsky,  but  the  resulting  relation  of  the  kernels  to  the  bicovariance  is 
essentially  the  same.  They  both  possess  a  finite  linear  relation  between  the  kernel 
function  and  the  bicovariance.  Accordingly  we  will  define  a  loose  class  of  models 
called  the  universal  bispectral  model  (UBM)  whose  members  have  this  characteristic. 
This  term  is  due  to  Wolinsky,  but  he  applied  it  to  only  one  specific  model.  The  main 
features  belonging  to  this  particular  model,  which  are  not  found  in  Wolinsky’s,  are 
that  its  kernel  function  can  be  made  to  be  symmetric  in  r  and  s,  i.e.,  g(r,s)  =  g(s,r); 
the  expression  relating  the  bicovariance  to  the  kernel  function  does  not  involve  any 
shifting  of  die  bicovariance;  and,  third,  the  power  spectrum  of  the  time  series  can  be 
changed  by  nxm  than  an  additive  constant  without  changing  the  bispectrum.  The 


ability  to  alter  the  power  spectnim  is  the  result  of  an  extra  degree  of  freedom  in  the 
kernel  function. 

From  the  noodel  equation  in  (1.1),  expressions  for  the  kernel  function  g(r,s) 
can  be  derived  in  terms  of  either  the  bispectrum  or  the  bicovariance.  These 
calculations,  which  are  done  in  the  case  of  discrete  time/continuous  frequency, 
constitute  the  heart  of  chapter  3.  In  chapter  4  the  subtleties  of  approximating  the 
discrete  tinte/continuous  fiequency  case  with  discrete  diiK/discrete  frequency  Fourier 
transforms  are  discussed  and  these  results  are  applied  to  some  specific  cases.  The 
discrete  time/discrete  fiequoicy  is  found  to  place  some  restrictions  on  the  arbitrariness 
of  the  bispectrum. 

There  are  two  applications  of  bispectral  inversion  that  the  writer  is  aware  of, 
but  potentially  there  are  many  more.  One  application  would  be  to  test  bispectral 
estimatcxs.  If  one  had  a  time  series  with  a  known  bispectrum,  one  could  then  test  an 
estimator  to  see  how  noise  affected  the  convergence  time,  the  amplitude,  or  the 
position  of  the  peaks  of  the  estimated  bispectrum.  Another  application  would  be  in 
the  area  of  communications.  A  communications  system  using  linear  signals,  for 
instance,  requires  that  multiple  transmissions  be  separated  in  either  a  time  window  or 
a  frequency  bandwidth.  In  other  words,  multiple  signals  of  distinct  frequency 
bandwidths  can  be  sent  at  the  same  time  and  still  be  recoverable,  and  multiple  signals 
occupying  the  same  frequency  bandwidth  can  be  sent  in  different  time  windows  and 
still  be  resolvable.  But  of  course  multiple  linear  signals  cannot  share  both  the  same 
frequency  band  and  the  same  time  window  and  still  transmit  recoverable  information. 
With  nonlinear  signals,  however,  this  is  not  necessarily  true.  Multiple  nonlinear 
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signals  occupying  identical  frequency  bands  and  occurring  simultaneously  can  be 
recovered  and  separated  provided  there  exists  a  nonlinear  characteristic  that  is 
different  The  bispectrum  is  such  a  characteristic.  Thus  by  exploiting  the  bispectral 
characteristics  of  the  signals,  a  communications  channel  could  simultaneously 
transmit  a  second  or  a  third  signal  in  the  same  frequency  band  as  the  first  and  still  be 
able  to  demultiplex  each  signal  without  amlsguity.  The  use  of  nonlinear  multiplexing 
to  encode  signals  is  the  topic  of  a  paper  by  Hinich  [1979].  However,  in  Hinich's 
paper  it  is  the  multiplexing' operation  that  is  nonlinear  and  not  the  signal  itself.  In 
either  case,  the  end  result  is  to  increase  the  number  of  signals  that  can  be  transmitted 
simultaneously  without  increasing  either  the  required  bandwidth  or  the  duration  of  the 
time  window  the  signals  occupy.  The  only  price  to  be  paid  is  the  additional 
ccxnputation  required  to  estimate  the  bespectium. 

Before  proceeding  to  the  derivation  of  the  g-kemels  it  would  be  well, 
considering  the  relative  newness  of  the  bispectrum  to  engineering  applications,  to 
present  a  review  of  the  bispectrum  and  its  characteristics. 


A 


CHAPTER2 


BISPECTRAL  TUTORIAL 


Essentially  the  bispectnim  is  one  component  in  the  set  of  statistical  quantities 
called  higher  order  spectra.  These  higher  order  spectra  are  the  frequency  domain 
equivalent  of  higher  order  statistical  quantities,  called  cumulants,  in  that  they  are 
Fourier  transform  pairs.  The  higher  order  spectra,  or  polyspectra,  provide  useful 
infoToadon  about  nonlinear  systems.  Just  as  linear  systems  can  be  analyzed  by  linear 
regression  techniques  using  the  power  spectrum,  so  in  an  entirely  analogous  manner 
nonlinear  systems  can  be  analyzed  by  multiple  regression  techniques  using 
polyspcctra.  In  the  event  there  is  only  one  time  series,  the  second  order  spectra 
reduces  to  the  ordinaty  power  spectrum,  while  the  third  order  polyspectra  is  called  the 
bispectrum.  If  there  are  two  time  series  the  above  quantities  become  the  cross 
spectrum  and  the  cross  bispectrum,  respectively.  Likewise  other  polyspectra  can  be 
defined,  such  as  the  fourth  order  spectra  which  is  commonly  called  the  trispectrum. 
However  the  higher  orders  are  increasingly  more  complex  and  difficult  to  actually 
compute.  Only  with  the  advent  of  faster  hardware  will  polyspectra  of  greater  order 
be  practical  to  compute.  In  addition  polyspectra  beyond  the  trispectrum  involve  more 
than  three  dimensions  making  them  veiy  difficult,  if  not  impossible,  to  visualize. 

The  purpose  of  this  chapter  is,  first,  to  provide  some  background  on 
polyspectra  culminating  in  a  definition  of  the  bispectrum,  and  secondly  to  intioduce 
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some  of  its  mathematical  properties,  paying  particular  attention  to  its  symmetries. 
The  impmtant  topic  concerning  estimation  of  the  bispectrum  will  be  passed  by 
without  much  comment  as  it  is  not  particularly  germane  to  this  thesis. 

2.1  Definitions:  Random  Processes,  Cumulants,  Moments  and 
Polyspectra 

Loosely  speaking,  detailed  knowledge  of  the  topics  in  measure  theory  on 
which  the  theory  of  random  variables  and  random  processes  rests,  and  from  which 
the  measure  theoretic  definitions  of  expected  value  and  other  statistical  quantities 
come,  is  not  necessary  in  the  actual  estimation  of  polyspectra.  However  it  is  good  to 
have  some  idea  where  these  concepts  come  from.  Hence  this  section  starts  with  a 
very  brief  definition  of  random  processes.  Further  details  can  be  found  in  works 
such  as  Priestly  [1981]  or,  for  the  more  ambitious,  Loeve[1977]. 

2.1.1  Random  Processes:  The  First  Building  Block 

A  random  process  is  defined  as  a  collection  of  Borel  measurable  functions 
(random  variables)  (  t  e  T  }  indexed  by  a  parameter  t  belonging  to  a  set  T,  and 

the  measure  is  some  probability  measure  defined  on  some  probability  space.  If  T  is 
the  set  of  real  numbers  {t:-«>o<t<oo),  denoted  by  SR,  then  the  process  is  called 
continuous.  If  T  is  the  set  of  natural  numbers  {  t :  ...  -2,-l,0,l,2. .  .  ),  denoted 
by  N,  then  (Xj)  is  a  discrete  parameter  process.  If  t  is  associated  with  time,  the 

process  is  called  a  time  series.  While  this  is  the  most  common  association,  it  is  not 
necessarily  the  only  one.  This  thesis  will  deal  primarily  with  a  discrete  parameter 
process.  It  is  worth  noting  that  there  are  two  parameters  involved  in  a  random 


process.  The  obvious  one  is  t,  and  die  not  so  obvious  one  is  the  subset  of  the  set 
Q  -  (ca :  OD  €  Q)  which  X  nu^is  Ui  Thus  fcH*  any  fixed  t  the  random  process 
reduces  to  a  random  variable,  while  for  any  fixed  (O  e  Q  the  random  process 
becomes  a  simple  function  of  t  sometimes  called  a  somp/e  parA  or  trajectory.  The 
reason  (O  is  not  commonly  written  will  be  clear  shortly. 

The  k  dimensional  probability  distribution  function  is  defined  as  a  function 
mapping  91^  to  St  such  that 

F  X(ti).  X(t2). . . .  X(tk)t*l““’^kl  =  X(ti)  ^  xj, . . .  ,X(ti.)  ^  xjt)  (2.1) 

where  P  is  the  probability  measure  mentioned  above  and  k  e  N  [Wise,  1986],  In  the 
foregoing  definitions  probability  spaces  and  measures  have  been  used,  but  have  not 
been  explicitly  defined.  The  reason  is  due  to  Kolmogorov  who  proved  that  given  a 
set  of  finite  dimensional  distribution  Unctions  satisfying  the  necessary  requirements, 
there  exists  a  probability  space  and  a  probability  measure  on  that  space  whose  family 
of  finite  dimensional  distribution  functions  coincides  with  the  given  family  of  finite 
dimensional  distributions  [Wise,  1986].  Thus  in  practice  it  is  not  necessary  to  show 
the  existence  of  probability  spaces  or  measures.  Having  the  distribution  function 
guarantees  their  existence.  Since  all  statistical  properties  are  completely  specified  by 
the  distribution  function,  neither  is  it  necessary  to  find  the  probability  measure  or 
space  in  order  to  rind  moments  or  cumulants.  It  is  for  this  same  reason  that  mis  not 
written  eiqplicitly. 

A  strictly  stationary  random  process  is  one  where,  for  k  e  N, 

F[  X(ti+t),  X(t2+t), . . .  ,X(tk+t)  1  =  F[  X(ti),  X(t2). . . .  X(tk)l.  (2.2) 


This  means  that  the  distribution  function  does  not  change  when  all  the  index 
parameters  t  are  shifted  an  equal  amount  up  or  down  the  axes  in  91^.  A  less  stringent 
restriction  is  stationarity  to  order  n.  A  random  process  is  considered  stationary  to 
\  (xder  n,  n  €  N,  if  Eq.  (2.2)  holds  for  any  k  ^  n. 

t 

2.1.2  Cumulants  and  Moments 


The  Hnal  definition  before  teaching  the  bispectrum  is  that  of  cumulants. 
Following  Rosenblatt  [1981]  we  define  the  moment  generating  function  (also  called 
the  characteristic  function)  of  the  vector  random  process  {Xaj(ti)...Xjil^(t|^)}  as 

k 

iTwX 

<!>(¥)=  <I>(¥p-Vfc)=Ele^‘  “)  (2-3) 


The  joint  moments  can  be  found  as  the  coefficients  of  (¥i>  •  ••¥  k)  i*'  the  Taylor  series 
expansion  of  <I>(y)  about  the  origin  provided  the  moments  are  finite.  As  Brillinger 
[197S]  points  out,  however,  this  is  not  really  a  problem  since  all  time  series  available 
for  processing  are  bounded.  The  cumulants  of  the  same  set  of  random  variables  are 
then  defined  as  the  coefficients  of  the  Taylor  series  expansion  of  the  log  d>(¥)  about 
zero.  Explicitly  this  gives  for  the  relation  of  the  joint  cumulant  sequence  of  order  r, 
where  r=ni+n2+  •  •  •  nj^  [NUdas  ef a/.,1987]: 


^  a{ln<I>(¥i,¥2-.¥k)} 

(ti,...tk)5(-j) - 

3¥i^2 


(2.4) 


v.=V2=-=yk=o 


The  tj  refers  to  the  random  variable  specified  by  evaluating  Xaj  at  tj.  Generally 
cumulants  (and  moments)  are  regarded  as  constants.  Random  processes,  however, 
are  a  collection  of  random  variables  of  which  each  combination  has  a  cumulant. 
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Hence  the  r***  order  cumulant  of  the  k-variate  random  process  can  be  viewed  as  a 
sequence  in  91^  indexed  by  the  same  k  parameters  as  die  random  vectm'  forming  the 
cumulant  or  moment.  The  corresponding  expression  for  the  joint  moment  sequence 
of  order  r,  which  is  defined  as 


i-  •  tk)  a  . .  x,^(tk)"^ ) 


(2.5) 


is  simply 


(<IK¥i.V2  -¥k)} 

(ti*tk)=  (-J)- 


a¥i  ^2- 


V,=yj=...=  y^0 


(2.6) 


where  aj.  j€  specifies  the  random  process  and  nj  specifies  the  power  to  which  the 
aj^  process  is  raised.  Further  examination  of  these  expressions  can  provide  some 
intuitive  explanation  of  the  difference  between  moments  and  cumulants.  Dr.  Choi's 
Ph.D.  dissertation  provides  good  insight  on  this  point  [Choi,  1984].  If  any  of  the 
random  variables  in  (2.3)  are  independent,  the  expectation  in  (2.3)  will  factor  into  the 
product  of  at  least  two  expectations.  Each  term  will  be  a  function  of  only  those 
random  variables  that  are  independent  of  all  the  other  random  variables.  This  is  to 
say  that  <I>(¥)  can  be  written  as  <lK¥)^i(¥i)^2(¥2)  'vhere  the  arguments  of  Oj  and 
O2  are  two  disjoint  sets  of  variables.  The  partial  derivative  in  (2.6)  is  not  affected  by 
this  factoring.  However,  when  the  logarithm  of  this  factored  expression  for  <IK¥)  is 
taken,  a  sum  will  result.  The  partial  derivative  in  (2.5),  to  which  the  cumulants  are 
related,  is  affected  by  the  factoring.  When  (2.3)  can  be  factored,  (2.5)  becomes 
identically  zero.  (  Notice  that  the  logarithm  does  not  commute  with  the  expected 
value  operator,  the  definition  of  a  cumulant  is  the  log  of  an  expected  value,  not  the 
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expected  value  of  a  log.)  Thus  while  the  cumulant  i»Dvides  information  about  the 
dependence  of  the  random  variables  [Choi,  1984],  the  moments  provide  informadon 
about  the  joint  distribution  functions.  This  information  is  the  statistical  analog  to  the 
information  provided  by  the  moment  of  inertia  about  the  distribution  of  mass  in  a 
physical  system.  Some  of  the  important  properties  of  cumulants  obtainable  from 
(2.3)  and  (2.4)  are  repeated  here  [Rosenblatt,  1981]: 

1)  C{Xi(ti),...,Xi^(^)}  is  symmetric  in  all  of  its  arguments. 

2)  If  any  group  of  X's  are  independent  of  the  others, 
C{X,(t,),...,Xk(tk))=0. 

3)  For  any  two  independent  groups  [Xj, ...,  X„)  and  [Xn+j,  ...,Xi[} 

C{Xi(ti)  +  X„+i(t„+i),....X„(t„)  +  X^(ty)]  =  C(Xi(ti) . X„(t„)) 

C(Xn^|(tu^j),...,X|[(t||r)). 

4)  For  a  set  of  constants  (ai...a)^}, 

C(aiX|(ti)...aitX|j(tjj))  *  ai...ai[C(Xi(ti)...X|j(ti^)). 


As  alluded  to  earlier,  the  polyspectra  are  defined  as  Fourier  transform  pairs 
with  the  cumulants.  However,  since  characteristic  functions  are  in  general  not 
available,  the  above  expressions  are  not  useful  in  the  actual  computation  of 
cumulants.  Since  they  are  also  not  directly  computable  from  representations  of  a 
random  process  [Kendall  et  al.,  1958],  it  is  of  interest  to  obtain  relations  for 
cumulants  in  terms  of  the  moments  which  are  directly  computable  from  ensemble 


representations  of  a  time  series.  These  relations  are  obtained  by  taking  the  logarithm 
of  the  Taylor  series  expansion  of  (2.3)  (  whose  coefficients  are  the  moments)  and 
expanding  it  again  as  a  Taylor  series  about  the  migin.  The  coefficients  of  this 
expansion  are  the  cumulants  by  definition.  This  is  worked  out  in  detail  in 
Kendall  etal.  [1958].  The  cumulants  can  then  be  found  by  suitable  matching  of 
terms.  An  expression  fcH’  the  moments  in  terms  of  the  cumulants  can  be  found  in  an 
analogous  manner.  This  derivation  was  done  for  the  general  case  by  Leonov  and 
Shriyaev  in  1959.  In  spectral  analysis  applications,  however,  only  the  simple 
moments  are  of  any  interest.  Simple  moments  are  those  where  n{sn2=...=n^=l. 
Hence  we  can  use  much  simpler  results,  which  are  only  vaUd  in  two  special  cases. 

Case  1: 

nj  s  n2  «...=  n|t  *  1  (simple  moments) 

<^ase2: 

n2  = ...  =  n^  =  0,  n|  arbitrary. 

With  these  restricdons  we  have  from  their  results; 

E{xi(t,). . .  Xk(tk)}=Xc(i)j)...c(t)^)  (2.7) 

where  \)  =  (vi...‘Ui(),  \>j  is  a  partition  of  the  integers  ( l,2,...,k},  and  the  summation 
is  over  all  such  partitions.  Likewise  the  expression  for  the  cumulants  in  terms  of  the 
moments  is 


c(xi(ti)...Xk(tk))=5^(-l)*^‘(p-l)!m(v,).„m(V 

u 


(2.8) 
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where  the  v  and  the  summatitHi  are  the  same  as  before  and  m  is  a  simple  moment 
about  the  origin.  Thus  simple  cumulants  of  order  k  can  be  expressed  as  polynomials 
of  simple  moments  of  mder  no  greater  than  k  and  the  moments  can  be  expressed  as 
polynomials  of  the  cumulants  of  order  no  greater  than  k.  The  information  contained 
in  the  excluded  cases  (Le.,  non-simple  moments)  is  not  being  thrown  away.  It  is 
actually  already  included  as  a  special  case  in  a  higher  order  (r^*'  order)  simple 
moment.  For  instance,  m{|2(ti.ti.t2)>  The  second  case,  i.e., 

(n2  =  ...  =  n^  =  0,  nj  arbitrary)  is  simply  a  special  case  of  the  first  one,  namely 
the  univariate.  In  addition,  for  a  stationary  time  series  it  turns  out  to  always  be  a 
constant  Hence,  without  much  danger  of  confusion  the  notation  can  be  simplified  by 
dropping  the  superscripts  since  they  are  always  one  in  the  cases  of  interest.  By  way 
of  illustration  consider  the  simple  fourth  order  cumulant  of  a  univariate  random 
process  (that  is,  X2=  ...  Xj.).  According  to  (2.8)  we  have 

Cllll(tl.t2>*3'*4)  =  ™llll(tl»*2>*3»*4)  ~  nill(ll.t2)™ll(l3»t4)  ~  niji(tj,t3)mii(t2,t4) 
-  mn(ti,t4)mii(t2,t3)  ~  m,(t,)miii(t2,t3,t4)  -  tn,(t2)mi  11(11,13,14) 
-mi(t3)mi  1  i(ti,t2,t4)  -  mi(t4)mi  1 1(1  i,t  2,13) 

+  2{mi(ti)mi(t2)mii(t3,t4)  +  mi(ti)mi(t3)mi  1(12,14) 

+  mi(ti)mi(t4)mii(t2,t3)  +  mi(t^mi(t3)mii(ti,t4) 


+  m,(t2)mi(t4)mi,(ti,t3)  +  mi(t3)mi(t4)mii(ti,t2)) 
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-  &ni(ti)mi(t2)mi(t3)mi(t4)  (2.9) 

If  we  now  assume  that  the  random  process  is  zero  mean,  i.e.,  m(t)=0,  and  that  it  is 
stationary  to  at  least  fourth  order,  (2.9)  becomes  much  simpler.  A  direct  consequence 
of  the  application  of  the  deHnidon  of  stadonarity  to  the  expected  value  operator 
(E[x]  =1  X  dF  )  is  that  moments  of  stadonary  random  processes  can  be  written  in 

terms  of  time  differences.  If  the  dme  lag  is  chosen  to  be  t„.i  =  t^-  tj  for  n 
€  (  1,2,.. .,k},  then  the  order  mcnnent  can  be  written  in  terms  of  k-1  lags 
because  the  lag  corresponding  to  n=l  has  become  idendcally  0  aiul  can  be  dropped. 
Again  the  notadon  can  be  simplified.  By  convention,  when  dealing  with  a  univariate 
sequence,  the  moment  we  have  written  as  m]],  j  is  often  written  as  m].^].^  . 

Following  this  convendon  (2.9)  simplifies  to 

-  m^x  -  m^x^m^Xy-^ ,)  -  m^x^jliX2-^  i)  .  (2. 10) 

Now  suppose  that  the  randtxn  variables  deHned  by  fixing  X2  and  X3  are  independent 
of  those  at  0  and  x^.  Then  the  expectation  of  the  first  term  on  the  right  side  of  (2.10) 
will  factor  making  m4(Xj,X2,X3)=m2(‘tj)m2(X3-^2)  *hc  first  two  terms  in  (2.10) 

would  cancel  out  The  last  two  terms  become  zero,  because  the  mean  is  zero,  causing 
the  whole  cumulant  to  vanish.  The  nonstadonary  cumulant  in  (2.9)  would  vanish  as 
well.  If  one  random  variable  were  independent  of  some,  but  not  all  of  the  others, 
then  some  but  not  all  the  terms  in  (2.9)  or  (2.10)  would  cancel.  Thus  the  cumulant 
measures  the  degree  of  dependence  between  the  random  variables.  Table  1  provides  a 
summary  of  the  reladon  between  cumulants  and  moments  up  to  order  four.  The  table 


is  written  in  the  univariate  notation,  but  can  be  generalized  to  the  multivariate  case  by 
simply  replacing  term  for  term  the  univariate  moment  (cumulant)  with  the  desired 
moment  (cumulant). 


Table  1  Relation  Between  (}umulants  and  Moments 

C|  =  mj  =0 


C2(ti)=  m2(Xi) 

Ci(Xi,Xj)=  m3(T,,t2) 

C4(Xi.T2  t3)  =  t04(Xi,X2X3)  -m2(Xi)m2(X3-X2) 

-m^Xj^mjUx^-x  i)-m2(X3)m2(X2-x ,) 


2.1.3  Polyspectra 

Using  the  results  of  the  last  two  sections  the  formal  definition  of  polyspectia 
follows  very  easily.  If  a  k-dimensional  random  process  satisfies  the  stationarity 
requirement  given  in  (2.2)  and  the  moment  niai...a]((h  *-V  exists,  then,  following 

Brillinger  [1965],  the  order  cumulant  can  be  written: 
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g^..  ...df »  (2.1 1) 


where  g,f..»^fi...fk)S(fi+ .  + f^)  is  the  ordtr  cumulant  spectral  density 

function.  Since  stationarity  causes  the  spectral  density  to  be  zero  except  when 
_  f  j-  {2- ...  -  ,  the  k*  order  density  can  be  written  as  a  function  of  k-1 

arguments  with  the  omitted  k***  argument  understood  to  be  the  sum  of  the  first  k-1 
frequencies.  The  expression  in  (2.1 1)  can  be  inverted  to  give  an  expression  for  the 
spectral  density  in  terms  of  the  cumulant;  or  equivalently,  the  k*  order  cumulant  is 
the  Fourier  transform  of  the  k***  order  polyspectra.  If  we  first  integrate  with  respect  to 
fj  using  the  dirac  function,  we  can  rewrite  (2.1 1)  as 


^  ^  •m 


(2.12) 


where  f^  is  not  written  but  is  understood  to  be  fk=  -  fj-  f2-  •••  -  fk-i  tutd  t„=tn-tk. 
Note  that  the  right  side  is  a  futKtitm  of  x  only.  This  produces  the  same  result  obtained 
earlier  for  stationary  time  series;  namely,  the  cumulant  can  be  written  in  term  of  time 
differences.  Thus  the  fiequency  domain  equivalent  of  writing  rrxrments  in  term  of 
time  lags  is  that  the  spectra  is  zero  except  at  points  where  all  the  fiequencies  sum  to 
zero.  Continuing  the  inversion  of  the  Fourier  transform  gives  for  the  spectral  density 
function  of  order  k-1  for  a  stationary  time  series: 


gtf  ••  ai^^l 


fk-i)-  •  jj  ' 


.(2.13) 
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Fw  the  univariate  case  with  k:s2  we  have  from  the  table  relating  cumulants 
and  moments  that  Cii(t)=m|i(T).  Since  for  a  real  time  series  m)i(x)  is  just  the 
autocorrelatimi,  we  see  that  the  second  order  cumulant  spectra  reduces  to  the  power 
spectrum  as  stated  earlier.  When  the  spectral  density  function  is  called  the 
bispectrum.  In  the  bivariate  case,  i.e.,  Cj2('f)  or  Cj  12(11, T2),  the  spectral  density  is 
called  the  cross-power  or  bispectrunL  In  this  thesis  we  will  be  interested  exclusively 
in  the  bispectrum.  In  particular  we  will  be  interested  in  the  discrete  frequency  / 
discrete  time  biqrectrunL  The  objectives  of  the  next  section  will  be  to  bridge  the  gap 
between  the  discrete  and  continuous  versions  of  the  bispectrum  and  to  present  some 
of  the  details  0[  the  bispectrum  that  will  be  important  in  the  inversion  process. 

2.2  The  Bispectrum 

In  concentrating  on  (2.13)  for  the  specific  case  where  k=3  we  will  make  some 
simplifying  assumptions.  Bom  this  point  onward  we  will  restrict  ourselves  to  real, 
stationary,  zero  mean  time  series.  The  assumption  of  reality  is  very  easily  justified. 
Any  experimentally  measured  time  series  is,  of  course,  real.  They  do  not  necessarily 
satisfy  the  other  two  assumptions  though.  The  zero  mean  assumption  is  easily 
achieved  by  simply  subtracting  the  mean.  The  stationarity  assumption,  on  the  other 
hand,  is  rarely,  if  ever,  completely  satisfied  in  experimental  data.  It  is  an 
approximation  that  is  valid  over  "short"  intervals.  But  in  the  case  being  investigated 
by  this  thesis  the  actual  model  generating  the  series  is  known  exactly.  Hence  we  can 
say  without  approximation  that  the  time  series  is  stationary.  With  these  assumptions 


let  us  first  rewrite  (2.13)  in  convenient,  more  common  notation.  Since  X(t)  is  real 
we  can  equate 
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C3(X|.T2)  =  =  E[  X(t)X(t-K,)X(t+t2]  /  R(Ti.T2) 


When  this  equivalence  is  applied  to  (2.13)  in  the  specific  case  of  the  bispectrum,  it 
gives: 


B(fi/2) 


=  jj  El  X(t)X(t-Ki)X(t4X2)]c*^^^’''^‘*^^dtidt2 


(2.14) 


where  B  is  the  bispectrum  and  the  integral  is  the  two-dimensional  Fourier  transform 
of  the  third  order  moment.  Fig.  2.1  shows  the  symmetry  regions  of  the  third  order 
moment  sequence,  C3(ti,T2)- 


8 


Equation  (2.14)  gives  a  pleasing  analytical  expression  for  the  bispectrum,  but 


it  is  not  the  nx>st  convenient  expression  to  use  to  compute  the  bispectrum.  A  better 
(Mie  can  be  obtained  by  going  back  to  (2.1 1)  and  inverting  the  k-dimensional  Fourier 
transform  directly  to  give 
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g3(fi.f3f3)5(fl  +  f2+f^ 


=  JJJ  C3(t^t3t3)  d^dtj  .  (2.15) 


Making  the  two  substitutions 

1)  C(t,.t2,t3)  =  E  [x(ti)X(t2)X(t3)] 

2)  x(t)  =  f  X(f)c^^’^df 


and  interchanging  the  order  of  integration  allows  (2.  IS)  to  be  rewritten  as: 

g3(fi.f2.f3)5(fi+f2  +  f3)  =  jjj  dfidfjdfjJJI'  dtidt2dt3 

X  E[x(f  ,)X(f2)X(f3)]  ^  (^2*  (frf2)h} 


(2.16) 


Carrying  out  the  the  three  integrations  with  respect  to  t  produces  the  product  of  three 

AAA 

Dirac  functions  6(fi-fi)6(f2-f2)5(f3-f3)  .  This  makes  the  integration  with 
respect  to  f  trivial,  giving  for  (2.16) 

g3(fi.f2.f3)5(fi  +  f2  +  f3)=  (2.17) 


Thus 


E[x(fpX(f^X(f3) 


=  0  unless  all  three  frequencies  sum  to  zero.  Taking  this 


into  account,  (2.17)  can  be  rewritten  to  give  the  more  revealing  expression  for  the 
bispectrum 


B(f,.f2)8f(0)  =  E 


X(f,)X(f2)X  (fi+f2) 


(2.18) 
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From  this  expression  one  can  see  that  the  bispectrum  is  nonzero  only  when 
the  fiequency  pair  (fj,  and  its  sum  are  coupled  so  as  not  to  be  statistically 
independent.  The  bispectrum  is  a  measure  of  quadratic  nonlinearity  because  it  is  a 
quadratic  interaction  between  f|  and  £2  which  will  produce  a  third  frequency  (fi+f2) 
that  is  statistically  dependent  to  the  first  two.  From  the  expression  in  (2.18)  it  is  a 
simple  matter  to  determine  the  symmetries  of  the  bispectrum,  although  they  can  also 
be  found  from  (2.14)  using  the  symmetry  of  the  third  order  moment  shown  in 
fig  2.1. 

2.2.1  Symmetries  of  the  Continuous  Frequency  Bispectrum 
From  (2. 1 8)  it  is  clear  that  f^  and  f2  are  interchangeable. 

Ptapcny  L 

B(fi,f2)  =B(f2,fi). 

This  produces  the  line  of  symmetry  labeled  1  in  fig  2.2.  Substituting  fj  (-fj) 
and  f2  -*  {-{2)  into  (2.18)  and  using  the  fact  that  X(0  =  X*(-  f)  gives  a  second 
symmetry. 

PcQpcnyZ: 

B(f„f2)  =B*(-fi,-f2)- 

This  produces  the  line  of  symmetry  labeled  2  in  fig  2.2.  Two  additional  relations  can 
be  found  by  first  substituting  in  fj  -♦  (fi+f2)  and  either  f2 (-  f2)  or 
f2  fl)  and  then  applying  the  conjugate  property  of  the  Fourier  transform. 


PtePfltyi; 


B(fif2)  =  B*(fi+f2.-f2) 


\ 


Application  of  property  2  to  each  of  the  terms  in  property  3  produces  two  more  terms. 
This  gives  five  terms  in  addition  to  property  1.  Applying  property  1  to  these  five 
gives  five  more  for  a  total  of  ten  terms.  These  symmetries  divide  the  bispectral 
domain  into  12  equivalent  regions  as  shown  in  fig  2.2.  Although  the  regions  do  not 
appear  to  be  of  equal  area,  each  region  does  contribute  equally  to  the  region  of 
support.  The  shaded  region  represents  the  region  of  support  for  the  bispectrum  of  a 
bandlimited  process.  The  bispectrum  is  uniquely  specified  by  knowing  the  values  in 
any  one  of  these  regions.  If  we  take  the  region  {(fi,f2):  fi  ^  f2;  f2^  0  }  to  be 
the  primary  region  (principal  domain)  then  the  1 1  other  equivalent  regions  are  given 
explicitly  in  terms  of  this  primary  region  by  (2.19).  Thus  for  a  bandlimited  process 
the  region  of  support  within  the  principal  domain  is  the  only  part  of  the  bispectrum 
that  needs  to  be  computed.  When  arbitrarily  specifying  the  bispectrum  it  should  be 
noted  that  the  bispectrum  must  be  purely  real  when  either  fi  or  f2  is  zero.  This  is 
evident  from  (2.18)  or  from  the  fact  that  this  is  the  only  way  to  satisfy  (2. 19)  on  the 
boundaries  between  regions  of  the  bispectrum  which  are  complex  conjugates  of  each 
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If  computers  worked  with  continuous  time  variables  we  could  stop  here. 
However  they  do  not.  As  will  be  seen  shordy,  this  discretization  produces  some 
significant  changes  in  the  the  frequency  plane  of  the  biqiectrum. 


BCfi.fj) 

□ 

B(f2.fi) 

m 

B  (-f2.fi+f2) 

B  (-fi,f|+f2) 

0 

B  (-frf2.fi) 

0 

B(-frf2.f2) 

0 

B*(-fi.-f2) 

0 

BVf2.-fl) 

0 

B  (f2.-f  rf2) 

0 

B(f,.-frf2) 

B*(fl+f2,-fl) 

0 

B*(fi+f2.-f2) 

12 

(2.19) 
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Fig.  2.2  Biapeetral  Symmetry  for  Continuoug  Frequency 


a)  Sfuided  area  represents  the  re0um  of  support  for  a  BandUmited 
Bispectrum  tuhose  highest  frequency  is 
B)  3aKfd  numBers  reference  esfpressions  in  'Ey.  (2.19) 
c)  The  mapping  of  an  orBitrary  magnitude  characteristic  is 
shoton  By  the  tu/elve  Boe(fs.  The  ’C‘  denotes  where  the 
Bispectrum  is  the  conjugate  of  the  principal  domain,  and  the 
triangfes  in  the  Bo)(fs  rtferenu  specific  comers  and  sides. 
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2.2.2  Symmetries  of  the  Discrete  Time/Continuous  Frequency 
Bispectrum 

The  discrete  tune  series  causes  the  bispectrum  to  be  periodic  in  fi  and  f2  with  a 
period  of  f,  or  1,  where  f,  is  the  sampling  frequency  and  1  is  in  terms  of  normalized 
frequency.  This  is  a  direct  consequence  of  the  periodicity  of  the  Fourier  transform  of 
a  discrete  time  series  as  can  be  seen  from  (2.18).  Thus  the  frequencies  that  had  to 
sum  to  zero  in  the  continuous  case  now  only  have  to  sum  to  0  (mod  1).  The  two- 
dimensitmal  FtMnier  transfemn  of  the  third  order  moment  (bispectrum)  now  becomes: 

B(fi.f2)=  2^2^  (2.20) 

—  ** 

where  the  domain  of  B  extends  over  the  whole  real  plane.  Although  B  is  defined 
over  the  whole  real  plane,  it  is  now  doubly  periodic  with  respect  to  both  f]  and  f2 
with  period  1.  This  implies  B(fi,f2)=  B(fi,f2+1)  =  B(fi  +  l,f2+l).  The 
fundamental  domain  is  defined  as  the  region  |(fl,f2)’  *^2^  ^  Yz’ 

the  principal  domain  continues  to  be  that  portion  of  the  fundamental  that  uniquely 
specifies  the  bispectrum.  The  tilda  is  to  distinguish  the  periodic  and  non-periodic 
versions  of  B.  Equation  (2.20)  is  the  discrete  equivalent  of  (2.14).  The  discrete 
equivalent  to  (2.15)  is  similar  except  diat  now  because  of  periodicity  (f]+  {2+  (3 ) 

may  sum  to  0(mod  1).  The  dirac  on  the  right  side  of  (2.15)  becomes 
5(n-[fi  +  f2  +  f3l)  where  n=  {0,±1,±2,®).  The  rest  of  the  development  of  the 

discrete  equivalent  to  (2.18)  follows  in  a  completely  analogous  manner  allowing 
(2.18)  to  be  written  in  the  discrete  case  as: 
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P 

f- 


B(f  1/2)  =  E[  X(fi)X(f2)X(n  -{f  1+  fj})]  (2.21) 

where  X(f)  =  5^  x(t)e 

tm.m 

where  the  overlined  portion  represents  the  modified  selection  rule  for  the  third 
frequency.  When  lie  inside  the  fundamental  domain,  the  only  two  values  of 

n  that  yield  a  third  frequency  that  is  also  between  -Y^  and  are  n  =  0,1.  The 

other  values  of  n,  which  put  the  third  frequency  outside  the  fundamental  domain, 
yield  only  redundant  infrmnation  since  X(ffn)  =  X(f).  When  n=0,  (2.21)  reduces  to 
(2.18):  the  resulting  symmetries  are  the  same  as  before.  When  n=l,  however,  some 
new  symmetries  result  Substituting  fj  -+  (1-  f|-  and  f2  -*  fi  or  {2  gives  two 
additional  pn^rerties. 

Additional  Symmetries  for  Discrete  Time: 

a) B(fi,f2)  =  B(l-fi-f2,f,) 

b) B(fi,f2)  =  B(l-fi-f2.f2)  . 

If  f]  and  f2  are  chosen  from  region  1  in  flg.  2.3,  these  additional  relations  simply 
repeat  the  previous  ones  in  (2.19).  If,  however,  they  are  chosen  from  region  2  there 
is  no  repetition;  they  yield  new  symmetries,  as  evidenced  by  the  fact  that  both  of  the 
resulting  terms  still  lie  within  region  1  on  fig  2.2.  It  now  appears  that  there  are 
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Fig.  2.3  Principal  Domain  of  PiBcrete  Time  Bispectrum 


additional  symmetries  within  region  1  on  fig.  2.2,  and  indeed  this  is  the  case  as  is 
shown  in  fig.  2.3a.  Applying  pr(^}erty  1  to  the  area  shown  in  fig.  2.3a  produces  the 
three  equivalent  regions  shown  in  fig.  2.3b.  The  use  of  property  2  will  give  six 
more  equivalent  regkms  in  the  fourth  quadrant  for  a  total  of  12  additional  symmetrical 
regions  as  shown  in  fig.  2.3b.  These  "extra"  regions  arise  because  of  the  double 
periodicity  oi  the  bispectrum.  Using  this  periodicity  these  12  "extra"  regions  can  be 
paired  with  the  original  12  regions  shown  in  Hg.  2.2  to  produce  a  different 


fundamental  domain.  This  alternate  fundamental  dcHnain  is  shown  in  fig.  2.4.  This 
is  the  fundamental  domain  that  would  result  from  the  application  of  (2.19)  to 
firequency  pairs  (fi^2)  in  both  regions  of  the  principal  domain  shown  in  fig.  2.3. 
This  is  the  more  usual  definition  oi  the  fundamental  domain.  It  is  simpler  in  that  it 
does  not  require  the  use  of  additional  symmetry  relations;  however,  the  two  are 
otherwise  entirely  equivalent  In  this  thesis  we  will  work  with  the  square  domain  for 
the  siiiq)le  reason  that  it  is  much  easier  to  perform  the  inverse  discrete  Fourier 
transfnm  over  a  square  region  than  over  a  hexagonal  one. 
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Notice  that  the  additional  symmetry  in  fig.  2.3  falls  outside  the  region  of 
suppoft  for  a  bandHmited  system.  This  means  that  for  a  pn^ioly  san:q>led,  stationary 
time  series  region  2  in  fig.  2.3  will  be  zero.  Recently  it  has  been  suggested  that  a 
nonzero  result  in  this  legitMi  inqrlies  either  a  nonstationary  or  aliased  time  series 
[Hinich  r/a/.,1988]. 

The  relations  defined  in  (2.18)  and  (2.14)  coupled  with  a  knowledge  of  the 
symmetries  given  in  section  2.2.2  provide  the  mathematical  framework  for  the 
inversion  of  the  bispectrum.  Much  detail  has  been  accorded  to  these  synunetries 
because  they  are  crucial  in  performing  the  inversion.  If  they  are  not  satisfied,  the 
third  order  nxxnent  obtained  by  transfimning  the  bispectrum  will  not  be  real  or  will 
not  have  the  prc^  symmetries  to  a  third  order  moment 
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Fig.  2.4 _ Fundamental  Domain  of  Discrete  Time  Bispectnun 


•Bo3(fid  numStrs  iiuCkate  equivalent  regions  of  Sisputrum 
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CHAFTER  3 


INVERSION  OF  THE  BISPECTRUM 


The  inversion  of  the  bispectnun  essentially  involves  finding  a  stochastic  time 
series  whose  bispectrum  satisfies  some  predetermined  constraints.  For  instance  it 
might  be  desired  to  produce  a  time  series  whose  bispectrum  is  a  constant  over  all 
points  or,  at  the  other  extreme,  the  requirement  might  be  to  have  a  bispectrum  that  is 
zero  at  every  point  but  one.  With  the  exception  of  two  philosophical  differences,  the 
method  used  to  invert  the  bispectrum  is  almost  the  exact  reverse  of  that  used  in  the 
more  usual  problem  of  estimating  the  bispectmm  of  a  time  series.  The  first  difference 
is  that  one  can  assume  to  start  with  the  true  bispectrum  and  not  just  an  estimate  of  it. 
When  attempting  to  estimate  the  bispectrum  from  a  finite  number  of  ensemble 
representations,  one  generally  assumes  to  have  only  part  of  the  time  series.  The  goal 
is  to  obtain  a  large  enough  subset  of  this  time  series  so  that  one  can  obtain  an  estimate 
of  the  bispectrum  to  within  some  acceptable  confidence  level.  In  the  present 
problem,  however,  we  assume  to  know  the  bispectrum  exactly.  But,  in  general,  the 
time  series  represented  by  this  bispectrum  (i.e.,  the  one  we  wish  to  compute)  is 
infinite.  In  a  finite  time  span  we  will  only  be  able  to  generate  a  finite  number  of 
ensemble  representations  of  this  underlying  time  series.  Hence  we  are  in  essence 
only  estimating  the  time  series.  The  quality  of  the  estimate  could  be  stated  in  terms  of 
the  number  of  ensembles  that  are  needed  to  produce  a  bispectral  estimate  that  is 
arbitrarily  close  to  the  true  bispectrum.  Since  this  number  could  vary  depending  on 
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the  estiiDator,  this  provides  a  means  to  perfonn  comparison  tests  between  bispectral 
estimators.  The  second  difference  is  that  while  every  time  series  has  only  one 
bispectrum,  there  are,  in  general,  many  different  time  series  that  all  have  the  same 
bispectrum.  This  is  completely  analogous  to  the  fact  that  different  time  series  (and  in 
general  an  infinity  d  time  series)  have  the  same  autocorrelation  function.  This  second 
difference  is  not  of  great  concern  here  since  all  that  is  needed  is  one  time  series.  The 
fact  that  there  may  be  mote  is  of  no  consequence. 

It  should  also  be  noted  that  the  results  of  this  chapter  are  valid  for  continuous 
time  without  any  further  modification.  The  problems  arise  when  one  attempts  to 
handle  discrete  frequencies,  as  will  be  seen  in  chapter  four. 

3.1  Necessary  and  Sufficient  Conditions  for  Invertibility 

To  invert  the  bispectrum  we  will  use  the  simplest  class  of  models  that  are 
ciqrable  of  reproducing  the  Inspectrum.  To  be  capable  of  reproducing  an  arbitrary 
bi^)ectrum  the  model  must  satisfy  three  general  requirements. 

1 )  It  must  have  a  non-zero  bicovariance. 

2)  The  resulting  expression  for  the  bicovariance  of  the  model  in  terms  of 
the  g-kemeis  must  be  a  function^  which  we  will  call  the  bicovariance 
function .  (The  domain  of  this  function  is  the  kernel  sequence,  and  the 
range  is  the  bicovariance  or  third  order  cumulant  sequence.) 

3)  The  domain  of  the  kernel  sequence  must  be  the  same  as  the  domain  of 
the  bicovariance. 
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The  first  one  is  obvious  and  included  only  for  completeness.  The  second 
condition  guarantees  that  the  kernels  will  be  mapped  to  only  one  bicovariance. 
Essentially  it  eliminates  the  possibility  that  two  different  bicovariances  will  both  have 
the  same  g  kernels,  and  ensures  that  it  will  be  possible  to  express  the  model  kernels 
in  terms  of  a  unique  bicovariance.  This  requirement  might  eliminate  some  AR 
processes.  Note,  however,  that  the  inverse  relation  need  not  be  a  function,  because  it 
is  permissible  to  have  different  kernel  sequences  produce  the  same  bicovariance;  the 
significance  of  this  will  be  apparent  later  when  we  attempt  to  invert  the  bicovariance 
function.  The  third  condition  guarantees  the  ability  to  reproduce  any  arbitrary 
bispectrum.  It  removes  any  constraints  on  the  bispectrum  or  bicovariance.  Without 
this  condition  it  would  be  possible  to  change  the  bicovariance  without  changing  the 
g  kernels.  There  are  many  models  that  could  be  constructed  which  fit  these  general 
requirements.  The  problem  is  to  find  a  model  which  has  a  bicovariance  function 
which  is  easily  invertible.  The  quadratically  nonlinear,  infinite  order  moving  average 
(MA)  model  of  the  universal  bispectral  class  is  ideally  suited  to  this  problem.  Since 
the  bispectrum  can  be  interpreted  as  a  measure  of  quadratic  coupling,  it  is  no  surprise 
that  the  model  best  suited  to  the  inversion  should  itself  be  quadratically  nonlinear. 
The  model,  introduced  in  chapter  one  as  the  universal  bispectrum  model,  is  repeated 
here  in  (3.1)  for  convenience. 

«•  eo 

x(0  =  n(0  +  +  r)Ttt  s)  (3.1) 

r=0x=0 

where  e(t)  and  Tj(t)  are  Gaussian  independent  identically  distributed  (i.i.d.)  random 
variables  with  zero  mean  and  variance  or  N(0,o^).  This  model  is  unique  in  that 
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the  tncovariance  function  is  a  finite  linear  relation  consisting  of  the  the  sum  of  six 
kernel  terms.  Models  outside  the  universal  bispectral  class  possess  more  complicated 
bicovariaiice  functicais  which  are  very  difficult,  if  not  impossible,  to  invert. 

3.2  Inversion  in  the  Time  Domain 

There  are  two  equivalent  methods  of  inverting  the  bispcctrum.  The  first  is  in 
the  time  domain.  The  three  main  steps  in  this  method  are  as  follows. 

1)  Solve  for  the  bicovariance  of  the  model  using  the  definition  of  the 
bicovariance. 

2)  Invert  the  resulting  expression  to  obtain  a  new  expression  for  the 
kernels  in  terms  of  the  bicovariancc. 

3)  Compute  the  g-kemels  from  the  relation  found  in  2)  thus  specifying 
the  time  series  by  means  of  (3.1). 

3.2.1  The  Bicovariance  Function  of  the  Universal  Bispectral 
Model 

The  first  step  in  obtaining  the  bicovariancc  is  to  substitute  (3.1)  into  the 
expression  for  the  third  order  cumulant  sequence  (or  the  moment  since  they  are 
equal): 

C3(Ti,T2)  =  E[X(r)X(r+T,)X(r+T2)]  .  (3.2) 


Doing  this  we  obtain: 
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C3l(T,,T2) »  E{€(t)  +  n(^)  +  ^  ^  girMt  +  r  )TJ(/  +  5) 

X  +rj)  +  Tflf  +rj)  +  ^ ^1+  +  ^1+ 


—  .  — 

X  e[f  +r2)  +  rKt  +^2)  ^  ^  ^  ^  ^ 


Expanding  this  product  yields  a  sum  of  27  tenns.  Invoking  the  fact  that  e(t)  and  ri(t) 
are  independent  and  zero  mean,  and  the  fact  that  odd  order  moments  of  Gaussian 
random  variables  are  zero,  reduces  this  expansion  to  only  six  non-zero  terms.  For 
example  the  first  term  of  the  expansion  would  be  E  |c(f )  £(t  +Tj)  €(r  . 

This  is  a  third  order  moment  of  a  Gaussian  random  variable  and  therefore  zero.  The 
term  containing  die  expression 


E 


€(/  +  r  )ij(/  +  j)£(r  +  T|+  r  Mt  +  Tj+  s)8(r  +  T2+  r  )T<t  +  T2+  s) 


is  also  zero  since  the  expectation  can  be  factored  into  the  product  of  two  terms  each 
of  which  is  a  third  order  moment  of  a  Gaussian  process.  Writing  only  the  six 
remaining  terms  gives  for  the  bicovariance 

PO, 

=  E[  £(r)  rff  +  t{)  ^  ^«(r^)£(x  +  T2+  r  )rj(r  +  T2+  s) 


«•  00 

+  £(f)  TK?  +  ^2)  ^  +  ^1+  +  ^1+ 
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+  rj(t)e(f  +  T^)  ^  +  ^2+  »")»<'  +  ^2+ 

—  — 

+  rj(t)e(t  +  T2)^  ■'■  ^*'*‘ '' 

>»  .  — 

+  tX^  +  ti)  fi(?  +  *2)  ^  +  r)iit+  s) 

—  <» 

+  £(f  +  T,)T)(r  +  Tj)  ^  Y^girMt  +r)tit+  i)] 


Taking  the  expected  value  of  this  expression  allows  us  to  write  (3.3)  as 


ciitpT^  *  <r^|  ^  ^girMs  +  V  t{)Sr  +  tj 
+  £  +  V  ^2>^'’  +  ^1) 

J—  — 

+  V  +  ■'^2) 

+  ^  -  Ti)^s  -  Tj) 


where  S(  )  is  the  Dirac  delta  function.  The  delta  function  collapses  the  double 
summation  to  just  a  single  term  allowing  one  to  write  (3.4)  as  the  linear  finite  relation 
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l'^2)  +  l)  +  l) 


+  g(-t‘^^  i-t-j)  +  r  V^2) 

This  is  a  very  simple,  pleasing  result,  but  it  can  be  simplified  even  more. 
Notice  that  this  expression  ( but  not  the  kernel  functitm  itself)  automatically  possesses 
the  symmetry  necessary  fcv  a  third  order  moment  of  a  real  time  series  since  the  terms 
on  the  right  side  are  exactly  the  same  as  the  symmetries  of  the  third  order  moment 
shown  in  fig  2.1.  Thus  substituting  into  (3.S)  the  symmetry  relations  for  a  third  order 
moment  will  produce  a  system  of  six  equations  each  containing  six  kernel  terms,  but 
in  each  case  the  right  side  will  remain  unchanged  as  expected.  But  the  symmetry  also 
has  the  less  desirable  effect  of  making  the  system  indeterminate;  another  method  of 
inversion  must  be  used.  One  simplification  comes  fiom  noting  that  for  certain  values 
of  T)  and  t2  not  all  six  g  terms  are  present.  Inside  any  of  the  three  regions  in  fig.  3. 1 

only  two  terms  from  the  right  side  of  (3.5)  exist.  On  the  boundaries  all  the  terms 
from  the  adjacent  regions  exist.  Thus  at  the  origin  {(r,s):  r=s=0)  all  six  terms 
contribute,  but  in  this  case  (3.3)  reduces  to  c^0,0)  =  6g(0,0),  On  the  other 

boundaries  only  four  terms  contribute.  Two  terms  are  excluded  because  the  values, 
for  which  the  delta  function  in  (3.2)  is  non-zero,  are  not  in  the  range  of  summation 
causing  the  whole  summatimi  to  vanish.  Exactly  which  terms  contribute  depends  on 
which  boundaries  are  being  considered.  Off  the  origin  and  boundaries  only  two 
terms  contribute  to  the  bicovariance.  The  other  four  do  not  contribute  for  the  reason 
given  above.  The  exact  relation  between  the  bicovariance  and  the  kernel  function  is 
given  in  (3.6). 
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Tj>0,  T2>0 

2g(T|,0)+g(0,Tj) 

Ti>0,  T2  =  0 

g(T,  -T2,'T2)  +  g(-T2.ti-r2) 

Ti>T2,  T2<0 

2g(-T,,0)+g(0,-ri) 

ri<0,  r2=Tj 

T2>Ti,  Ti<0 

2g(T2,0)+g(0,T2) 

T2>0,  ri  =  0 

6g(0,0) 

o 

II 

o 

II 

(3.6) 


This  simplification  is  the  result  of  causality  and  therefore  could  also  have  been 
obtained  from  (3.S)  directly  by  noting  that  g(r,s)  =  0  when  r  <  0,  or  s  <  0. 

3.2.2  Inversion  of  the  Bicovariance  Function 

The  next  step  in  the  inversion  of  the  bispectrum  is  to  invert  (3.6)  and  express 
the  kernels  in  terms  of  the  bicovariance.  However  the  inverse  of  this  equation  is  not  a 
function.  There  are  an  infmity  of  kernel  functions  that  will  satisfy  (3.6).  This  arises 
from  the  fact  that  in  order  to  uniquely  specify  the  bicovariance  we  only  need  to  know 
the  values  of  C3(r,s)  in  the  region  ((r,s):  s  >  0,  r  >  s)  ).  For  any  point  in  this 
region,  however,  there  are  two  kernel  terms  that  contribute  as  can  be  seen  in  fig.  3.1. 
This  amounts  to  an  extra  degree  of  freedom.  The  exua  degree  of  freedom  is  actually 
not  a  problem,  but  rather  it  is  an  advantage.  We  are  free  to  specify  any  relation 
between  g(r,s)  and  g(s,r)  that  we  desire.  (Note  that  we  are  not  talking  about  the 
specific  term  g(r,s),  but  rather  how  the  weights  in  the  region  {(r,s):  s  >  0,  r  >  s) 
are  related  to  the  weights  in  the  region  {(r,s):  r  >  0,  s  >  r)  as  shown  in  fig.  3.2  ). 
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Fip.  3.2  Relation  Between  the  Kernels 

Q(r.s) 


Since  the  sum  cf  ^p,q)  and  g(q,p)  cmtriBuUs  to  c(p^), 
the  vahus  of  g(p4)  aitdg(q,p)  an  arhitrary  as  tong  as  the 
sum  is  the  same. 


The  bicovariance  in  any  legicm  is  dependant  (mi  only  the  sum  of  two  kernel 
terms.  This  nneans  the  value  of  the  individual  terms  can  vary  without  affecting  the 
bicovariance  as  long  as  the  sum  remains  constant.  The  extra  degree  of  freedom  is 
one  reason  this  particular  model  was  used  instead  of  the  one  suggested  by  Wolinsky 
[1988].  These  two  degrees  of  freedmn  give  one  the  ability  to  generate  multiple  time 
series  which  have  identical  bicovariances,  but  whose  other  moments  (and  hence 
polyspectra)  ate  in  general  different.  This  is  because  the  other  moments  will  in 
general  be  dependent  on  the  value  of  individual  kernel  terms  and  not  just  on  the  sum 
of  two  terms.  This  will  be  explored  Iniefly  when  we  look  at  the  autocorrelation  (rf*  the 
model.  Fot  the  purpose  of  reproducing  the  bispectrum,  let  g(r,s)  =  g(s,r).  This 
provides  the  simplest  possible  relation  f(V  (3.6).  We  may  now  write: 
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g(^l.t2)  = 


ti>0,  X2>0 

ti>0,  X2=0 

^C3(0.0) 

tl  =  0,  t2=0 

(3.7) 

The  rest  of  (3.6)  does  not  need  to  be  inverted.  The  other  expressions  only  duplicate 
the  information  given  by  (3.7).  This  should  be  no  surprise  since  the  weights  are  zero 
outside  of  the  domain  defined  in  (3.7). 

Given  a  bicovariance  one  could  at  this  point  compute  the  kernels  from  (3.7) 
and  then  compute  a  time  series  using  (3.1)  that  had  the  given  bicovariance.  As  it 
turns  out,  this  is  a  rather  slow  way  to  compute  the  time  series.  It  is  much  faster  to 
compute  the  Fourier  transform  of  the  time  series  and  then  take  the  inverse  transform. 
Since  the  goal  is  to  compute  a  time  series  given  its  bispectrum,  it  would  also  be  nice 
to  stay  entirely  in  the  frequency  domain  and  be  able  to  go  directly  from  the  bispectrum 
to  the  Fourier  transform  of  the  time  series  without  the  need  to  inverse  transform 
continuous  frequency  quantities  into  the  time  domain.  Having  to  relate  the  model 
parameters  to  the  bicovariance  in  the  time  domain  causes  approximations  to  enter  the 
picture.  The  approximations  arise  from  representing  the  bicovariance,  which  is  the 
inverse  transform  of  a  continuous  function,  by  the  inverse  discrete  Fourier  transform. 
Thus  staying  in  the  frequency  domain  would  in  effect  allow  one  to  know  the  Fourier 
transform  of  the  time  series  exactly.  Thus  one  would  only  have  to  resort  to  an 
approximation  when  computing  the  time  series  from  its  transform.  But  it  turns  out 
that  this  trouble  spot  is  not  evaded  so  easily.  The  reason  is  that  one  cannot  exploit 
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the  causality  in  die  finqueocy  domain.  Hius  there  is  no  waytoinvert  wliat 

would  be  the  Fourier  tiansfonn  o£  the  hicovariance  function  except  in  the  non-causal 
case  where  the  kernels  equal  the  tticovariance  scaled  by  a  factOT  of  six.  However,  it 
is  still  informative  to  look  at  the  inversion  process  in  the  frequency  domain  because, 
having  found  the  kernels  in  the  time  domain,  they  can  then  be  transformed  to  the 
frequency  domain  where  the  computation  of  the  time  series  is  more  easily 
accooqilished.  But  before  doing  this  we  will  look  briefly  at  the  autocorrelation  of  the 
model. 

3.2.3  Auto>correlation  of  the  Universal  Bispectral  Model 

The  procedure  for  conqiuting  the  autocmrelation  of  the  model  is  entirely 
analogous  to  what  was  done  eariier.  Using  the  definition  C2(t)  n  E  [  xCOxCt-K}] 
gives  for  the  second  order  cumulant  sequence: 

— 

C2!(t)  =  2c^fi[t)  +  O^X  S  (3  8) 

w=MB(rn) 

It  is  not  obvious  that  (3.8)  has  the  necessary  symmetry,  but  it  does.  The  important 
thing  is  that  this  equation  is  sensitive  to  the  relation  between  g(r,s)  and  gCs^*).  Thus, 
as  claimed,  it  is  possible  to  produce  different  autocorrelation  functions  without 
changing  the  bicovariance  simply  1^  changing  the  relation  between  g(r,s)  and  g(s,r). 
However,  it  is  certainly  not  a  trivial  problem  to  determine  the  relation  necessary  to 
produce  a  desired  autocotrelaticm  function.  For  cmnpleteness  the  power  spectrum  is 
given  below. 
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G(fi,  f2  )  is  the  two-dimensional  Fourier  transfonn  of  the  kernels.  The  limitation  on 
arbitrarily  specifying  the  power  spectrum  stems  from  the  fact  that  once  the  bispectrum 
is  picked  to  be  nonzero  at  a  certain  point,  say  (fi/2)>  the  power  spectrum  cannot  be 
set  to  zero  at  fi,  fa  or  fi+  fa.  In  the  non-causal  model  the  transfonn  of  the  weights  is 
simply  the  bispectrum  scaled  by  a  factor  of  six.  This  will  be  shown  in  the  next 
section.  But,  as  will  be  shown,  even  in  the  causal  case  the  modulus  of  G(fi,  fa) 
bears  some  resemblance  to  the  modulus  of  the  bispectrum.  Thus  it  is  an  interesting 
sidelight  to  note  how  the  power  spectrum  is  related  to  the  bispectrum  in  the 
non-causal  case.  This  is  shown  in  fig.  3.3.  In  the  figure  the  diagonal  lines  represent 
the  path  of  the  integral  for  specific  values  of  f.  The  bispectral  synunetiy  forces  the 
power  spectrum  to  be  symmetric  about  the  origin  as  it  should  be. 
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3.3  Inversion  in  the  Frequency  Domain 

The  first  item  on  the  agenda  is  to  find  the  Fourier  transform  of  the  model.  We 
first  substitute  into  (3.1)  the  Fourier  integral  for  e(t+r)  and  T|(t+r).  This  gives  for 
x(t) 

x(/)  =  7K0  +  e(/)+2&'’’^>l  •  (3.9) 

faOjaO  J  l/  J_\/ 
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Interchanging  die  order  of  integradon  and  summation  and  leananging  the  terms  gives 


x(t)  =  rK0  +  £(0  +  JJ 

~/i 


fl^xT  +/2*) 


'(3.10) 


Since  the  kernels  are  zero  for  r  ot  s  <  0,  the  summation  on  the  right  constitutes  the 
conjugate  of  the  Fourier  transform  of  g(r,s)  or  G*(fi,fi).  Taking  the  Fourier 
transform  of  both  sides  of  (3.10),  and  bringing  the  Fourier  transform  operator 
throu^  the  integrals  gives 


X(/)  =  tK/)  +  £(/)  +  Jj*' 


(3.11) 


where  <j(f  1/2)=  Y, 


•flmfir  +/2J) 


r  9e  ^ 


The  transform  of  the  exponential  is  simply  the  delta  function.  The  delta  function 
allows  evaluation  of  one  integral,  reducing  the  double  integral  to  a  single  integral  and 
allowing  us  to  write  (3.1 1)  as 


X(/)  =  »?(/)  +  £(/)+  I  <ri^i)tl(f-/i)GVi/-/i)  . 


(3.12) 


At  this  point  we  can  substitute  (3.12)  into  (2.18).  This  gives: 


mx4i>  =  £  I  (n^i)  +  ^i)  +  ^vfSx  -/)  C  V./i  -/)) 

f‘A 

X  {Tj(/-2)+fi(ir2)+  (r0(nw2-f)G*(f,f2-f)) 

hA 

X  (Tj(fi+/2)  +  ^l+/2)+  4^£W»7ai+/2-/)cV./i+/2-/)}l  •  (3.13) 

/  •/ 


What  follows  now  is  similar  to  what  was  done  in  the  time  domain  except  that  here  we 
are  computing  the  bispectrum  directly  in  terms  of  the  two-dimensional  Fourier 
transfonn  of  the  kernels  instead  of  the  bicovariance  in  toms  of  the  kernels.  There  are 
27  terms  in  the  expansion,  and  again  only  six  remain  after  taking  the  expected  value. 

Note  that  ^£(jri)fl(^2)]  =  O  5(^i+/2)  allows  us  to  write  for  the  bispectrum 


Bffx/2) 


=  a^['V[c?(r./i+/2- 
Jv 


+  G  <f,frf)S(f-*-fx)  +  G  {f,fx-f)S(f^fi) 


(3.14) 


The  Dirac  functions  allow  us  to  evaluate  the  integral;  thus  we  can  write  the  bispectrum 
in  terms  of  the  transform  of  the  kernels  as 


B(fl/2)=<^  iG(fl/2)  +  GCr2/l)  +  G  (-/2./iif2)+G 
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+  G*(/i¥>*/a>+  i) }  .  (3. 15) 

This  relation  is  veiy  similar  to  (3.5)  in  form,  and  in  fact  is  simply  the  term  by 
term  two-dimensional  Fourier  transform  of  (3.5).  This  expression  cannot  be 
inverted  as  easily  as  (3.5)  was,  because  the  causality  of  the  kernels  cannot  be 
exploited  as  it  was  in  the  time  domain  expression.  However,  using  the  knowledge 
gained  by  the  inversion  in  the  time  domain,  one  can  obtain  an  expression  for 
G(fi,f2)  in  terms  of  B(fi,f2).  While  this  result  will  not  be  very  useful  from  a 
computational  aspect,  it  is  helpful  in  providing  an  intuitive  understanding  about  the 
lelatitm  between  the  kernels  and  the  bispectrum.  We  can  rewrite  (3.7)  as 

g(rj)  =  j  c(rMr)C/(s)  +  c(r^)S(s)C/(r)  +  j  c(r,s)C/(s)t/(r) 

-|c(r,j)5(s)fi(r)  (3.16) 


where  5(.)  is  the  Kronecker  delta  function  and  U(.)  is  the  discrete  equivalent  of  the 
Heavyside  step  function.  Taking  the  two-dimensional  Fourier  transform  of  both 
sides  gives  G(fi,f2)  on  the  left  side  and  some  function  of  the  bispectrum  on  the  right 
side.  The  step  and  delta  functions  by  which  the  bicovariance  is  multiplied  can  be 
regarded  as  separable  two-dimensional  functions.  Thus  their  two-dimensional 
Fourier  transform  is  simply  the  product  of  each  of  the  one-dimensional  Fourier 
transforms.  G(fi42)can  now  be  easily  written  as; 


48 


J'2jfi 


-  |Wl./2)*l,,  ,7^ 
Mi  M2  (3-17) 


(2«)/Lf2 


where  the  asterisk  denotes  a  two-dimensional  convolution.  One  can  see  why  this 
would  be  a  difricult  relation  to  use  for  computation,  but  it  does  tell  us  that  the 
transform  of  the  kernels  are  a  "smeared"  version  of  the  bispectrum.  They  bear  some 
resemblance  to  the  bispectrum.  It  is  worth  mentioning  that  in  the  non-causal  case 
alluded  to  earlier,  the  transform  of  the  kernels  is  simply  the  bispectrum  scaled  by  a 
constant. 

This  chapter  has  described  the  mathematics  involved  in  the  inversion  of  the 
bispectrum  in  the  case  of  discrete  time  domain  quantities  that  are  inrinite  in  duration 
and  frequency  quantities  that  are  continuous  in  nature.  Of  course,  in  order  to  actually 
implement  the  methods  set  forth  in  this  chapter  one  must  deal  with  finite  time 
sequences  and  discrete  frequency  functions.  This  means  that  the  bispectrum  will  be 
approximated  by  samples  from  a  finite  number  of  points,  and  the  bicovariance  will 
be  found  by  applying  the  inverse  discrete  Fourier  transform  to  this  approximation. 
The  implications  and  restrictions  inherent  in  this  approximation,  as  well  as  ways  to 
overcome  them,  will  be  the  subject  of  the  next  chapter. 


IMPLEMENTATION:  METHODS  AND  RESULTS 


The  purpose  of  this  chapter  is  twofold.  The  Hrst  is  to  discuss  the  subtleties 
involved  in  the  implementation  of  the  inversion  methods  derived  in  the  previous 
chapter.  These  subdeties  revolve  around  the  fact  that  we  are  discretizing  a  continuous 
quantity  (frequency)  in  order  to  compute  transforms  via  the  DPT.  While  to  some 
such  things  are  mere  trivialities,  they  nevertheless  are  necessary  if  one  is  to  obtain 
correct  results.  The  second  purpose  for  this  chapter  is  to  give  the  step  by  step  results 
of  an  inversion  of  the  bispectrum. 

4.1  Obtaining  an  Unaliased  Bicovariance 

Up  to  this  point  we  have  assumed  the  bispectrum  to  be  a  continuous,  periodic 
function.  Its  inverse  transform  is  found  by  an  integration  over  the  fundamental 
domain.  We  now  want  to  sample  this  continuous  function  and  use  the  inverse 
discrete  Fourier  transform  (DFT)  to  compute  the  bicovariance.  Since  we  presume  to 
know  the  continuous  Inspectrum  we  are  free  to  sample  with  any  frequency  resolution 
we  desire.  The  resolution  that  is  used  will  place  a  constraint  on  the  minimum  length 
of  the  sequence,  but  otherwise  does  not  present  any  fundamental  problems.  The 
fundamental  difficulty  is  in  the  sampling  itself.  It  causes  the  bicovariance  to  be 
represented  as  a  finite  two-dimensional  sequence.  This  flniteness  raises  an  immediate 
problem.  Since  the  bicovariance  is  in  general  defined  in  two  dimensions  over  the 
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interval  (-  «>, «»),  the  finite  sequence  obtained  firtmi  the  inverse  DFT  of  the  sampled 
bispectnun  is  at  best  an  aliased  version  of  the  bicovariance  and  not  the  quantity  with 
which  the  kernels  are  defined.  definiteness,  suppose  the  bispectrum  is  sampled 
to  form  an  [2N-1  x2N'l]  sequence,  thereby  constraining  the  bicovariance  to  be  an 
[2N-1  X  2N-1]  sequence  defined  over  lags  in  the  interval  [-N+1,  N-l],  For  any 
arbitrarily  specified  bispectrum  Uiere  is  no  a  priori  method  to  determine  in  advance  if 
the  bicovariance  at  lags  greato’  than  ±(N-1)  is  relatively  small  or  not.  If  it  is  small, 
the  error  due  to  aliasing  may  be  negligible.  But  generally  one  cannot  assume  this. 
The  only  way  to  determine  aliasing  is  to  simply  compute  the  inverse  transform  of  the 
bispectrum.  By  looking  at  the  bicovariance  outside  the  region  of  support  one  can 
determine  whether  it  is  aliased  or  not.  The  values  outside  the  region  of  support  for  a 
finite  record  but  within  the  [2N-1  x  2N-1]  square  will  be  zero  if  the  bicovariance  is 
unaliased.  These  values  are  represented  by  the  empty  squares  in  fig.  4.1.  This  is 
because  the  bicovaiiance  in  these  regions  is  equivalent  to  the  bicovariance  in  regions 
where  at  least  one  of  the  lags  is  greater  than  ±(N-1).  Thus  a  nonzero  value  in  this 
region  inqrlies  diat  the  bicovariance  extends  farther  than  the  space  allotted  to  it  by  the 
DFT.  One  can  easily  verify  this  by  choosing  a  point  in  the  region  containing  empty 
squares  in  fig.  4.1  and  applying  die  syrmnetry  relations  of  the  bicovariance  to  it. 
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Pier.  4.1  Zero  and  Non-Zero  Reerions  of  a 

Finite  Bicovariance 

C(t,s) 

K-i 

i  s 

a  a  a  a  a  a  a 

■  ■■■■■  ■ 

a  a  a  a  a  a  m 

■  ■■■■■■ 

a  a  a  a  a  m  m 

■  ■■■■■  ■ 

a  a  a  a  m  m  m 

■■■■■■■ 

a  a  a  m  m  m  m 

■  ■■■■■  ■ 

■  ■■■■■■ 

■  ■■■■■■  r 

■  ■■■■■ 

■  ■■■■□  □ 

■  ■  ■  ■  a  D  □ 

■  ■  ■  □  a  □  a 

■  ■  a  a  a  □  □ 

■  a  a  □  a  □  □ 

o  o  a  D  a  a  □ 

■9Ul 

.  T 

■  . . .  iPoint  wfitft  Biccvariaux  is  nonzero 

□  . . .  Taint  roture  unaGastd  bicovarianu  is  zero; 

but  an  aGased  bkovariance  is  nonzero 

When  estimating  in  the  forward  direction,  i.e.,  time  series  to  spectra,  a  time 
domain  window  is  used  to  truncate  the  bicovariance  in  a  non-abrupt  fashion.  The 
assumption  is  that  the  bicovaiiance  is  negligible  for  lags  beyond  a  certain  value.  This 
is  essentially  what  is  needed  here.  But,  while  applying  such  a  window  is  a  relatively 
simple  process  in  the  time  domain,  it  is  quite  the  opposite  i.i  the  frequency  domain. 
There  it  involves  a  convolution  with  the  transform  of  the  window.  But  that  is  the 
least  of  the  problems.  A  bigger  concern  is  that  one  does  not  know  the  duration  or 
shape  of  the  bicovariance.  Thus  it  is  impossible  to  apply  any  sensible  criteria  in 
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picking  a  suitable  truncation  boundary.  One  caimot  determine  die  lags  beyond  which 
the  tncovariance  is  negligible.  The  result  is  that  although  a  window  is  needed,  one 
does  not  have  any  way  of  knowing  or  sensibly  choosing  what  window  parameters  to 
use;  specifically,  how  long  shmild  the  window  be  to  include  all  the  significant  points 
in  the  bicovatiance.  However,  there  is  another  way  to  look  at  this  problem.  Let  us 
leave  the  windowing  approach  for  the  moment  and  pursue  the  problem  from  a 
different  perspective. 


Another  method  of  obtaining  the  unaliased  bicovariance  is  to  find  a 
bispectrum  which  has  a  finite  bicovariance.  This  is  not  as  hopeless  as  it  may  initially 
appear.  Since  the  bicovatiance  can  be  arbitrarily  defined  over  the  interval  (-<»,  oo), 
there  is  no  reason  why  it  could  not  he  zero  for  lags  greater  than  some  L. 


at  1.^2)  = 


0 


-L  ^  Ti,T2  ^  L 
elsewhere 


(4.1) 


This  places  a  restriction  on  the  arbitrariness  of  the  bispectra  that  can  be  reproduced. 
One  can  now  consider  only  those  bispectra  which  possess  finite  bicovariances.  This 
may  seem  Uke  a  very  severe  restriction,  but  this  exact  assumption  is  commonly  made 
when  estimating  the  bispectrum  from  a  time  series.  Nevertheless  it  is  the  price  that 
must  be  paid  for  the  use  of  the  DFT.  Fortunately  this  does  not  place  severe  practical 
limitations  on  what  sorts  of  bispectra  can  be  reproduced.  It  turns  out  that  the  bigger 
problem  is  not  only  finding  such  a  bispectra,  but  also  finding  one  whose  magnitude 
characteristics  are  close  to  those  of  the  desired  bispectrum. 
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A  word  should  be  inserted  here  to  avoid  possible  confusion  when  attempting 
to  apply  what  has  been  said  i^re  to  the  analogous  case  of  the  power  spectrum. 
Unless  one  maintains  a  clear  distinction  between  the  estimate  of  a  quantity  and  that 
quantity  itself,  ctmfusion  may  easily  result.  For  instance,  it  is  well  known  that  the 
inverse  Fourier  transform  of  the  power  spectrum  produces  an  aliased  autoconelation 
function.  This  statement  is  referring  to  the  situation  where  one  has  an  M  point 
random  time  series  belonging  to  a  random  process  and  possessing  a  coherence  time 
much  less  than  M  (or  at  least  half  of  M).  Thus  one  assumes  that  the  true  correlation 
only  possesses  significant  values  for  time  lags  in  the  interval  {t:  -M/2  ^  x  ^  M/2} 
thereby  allowing  the  power  spectrum  to  be  estimated  from  an  M  point  DFT. 
However  the  estimate  of  the  autocorrelation  sequence  obtained  from  one  M  point 
ensemble  is  2M-1  points  long.  This  covariance  estimate,  unlike  the  true  covariance, 
is  of  course  not  obtainable  from  an  M  point  estimate  of  the  power  spectrum.  On  the 
other  hand,  if  one  assumes  to  know  the  true  power  spectrum  representing  the  true, 
finite  M  point  autocorrelation  function,  then  one  can  obtain  it  from  the  M  point 
inverse  transform  of  the  powo*  spectrum.  In  the  problem  that  is  the  subject  of  this 
thesis  we  assume  to  know  the  true  bispectra;  the  problem  is  finding  one  which 
actually  has  a  finite  bicovariance. 

Finding  a  bispectrum  which  has  a  fmite  bicovariance  is  conceptually  simple, 
but  not  as  easily  implemented.  In  cedar  to  produce  a  finite  bicovariance,  one  merely 
has  to  satisfy  certain  additional  constraints  in  the  bispcctrum.  These  additional 
relations  are  found  by  simply  setting  the  DFT  expression  for  a  pardcular  value  of  the 
bicovariance  equal  to  zero.  For  example,  if  N=:32  then  C(20,-20)  is  outside  the 
non-zero  region  of  an  unaliased  Incovariance.  Setting  this  equal  to  zero  gives 
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/i-o  A-o 


In  general  there  are  1/4(N2-1)  such  points  in  the  bicovariance.  This  gives  the 
following  linear  system  of  equations: 


A-o  A=o 

where  ie  {1,. .  .1/4(A^-1)} 


(4.2) 


Even  for  small  N,  solving  this  system  of  equations  is  no  small  task.  It  is 
further  complicated  by  the  fact  that  die  system  is  over  determined,  forcing  one  to  pick 
the  other  3/4(N2-1)  bispectral  points  before  solving  the  system.  These  have  to  be 
chosen  in  such  a  way  that  the  resulting  bispectrum  has  magnitude  characteristics 
similar  to  the  desired  bispectnim.  Since  the  equations  bear  close  resemblance  to  the 
DPT  and  in  fact  differ  only  by  a  constant,  one  can  solve  them  numerically  by  an 
iterative  method  using  the  DFT.  This  is  accomplished  by  specifying  the  bispectrum 
over  every  other  point  in  a  two-dimensional  fashion,  taking  the  inverse  DFT,  and 
then  setting  to  zero  the  proper  points  in  the  bicovariance.  The  modified  bicovariance 
is  then  transformed  back  to  the  bispectrum  where  the  initially  specified  points  are 
reset  and  the  other  points  are  left  alone.  This  process  is  continued  until  the  desired 
degree  of  accuracy  is  reached.  The  iteradon  will  converge  as  long  as  there  are  at  least 
1/4(N2  - 1)  points  left  unspecified  in  the  bispectrum.  In  other  words,  all  diat  is  being 
done  is  to  trade  the  choice  of  setting  values  in  the  bispectnim  fm-  the  choice  of  setting 
an  equal  number  of  values  in  the  bicovariance.  In  this  way  one  can  modify  a 
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biq)ectrum  to  force  it  to  possess  a  finite  bicovariance.  By  ^reading  the  specified 
points  out  in  the  fashion  indicated  above  we  are  enaUed  to  satisfy  both  the  symmetiy 


lequiiements  of  the  bispectnun  and  the  goal  diat  the  nxxlulus  of  the  modified 
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bUpectnun  have  close  agreement  with  the  desired  bispectnun.  The  end  result  of  this 
iteration  is  singly  the  equivalent  ttf  windowing  in  the  time  domain,  one  does  not 
window,  then  there  will  be  nonzero  values  in  the  part  of  the  bicovariance  shown  by 
I  empty  squares  in  fig.  4.1.  Blindly  proceeding  with  the  inversion  is  equivalent  to 

>  setting  these  values  to  zero.  This  course  changes  the  biq)ectnim  in  an  unknown 

and  hence  uncontrollable  way.  Use  of  the  iteration  algorithm  also  changes  the 
bispectrum,  but  it  does  so  in  a  way  diat  allows  one  to  control  what  gets  changed  and 
what  does  not.  Finally,  it  is  worth  noting  that  the  surface  integral  of  the  modified 
bispectrum  over  the  fundamental  domain  is  unchanged.  In  other  words  the  points 
added  to  the  modified  bispectrum  contribute  in  such  a  way  so  as  not  to  change  the 
skewness  of  the  time  series. 

4.2  The  DFT  of  the  Model 

The  purpose  of  this  section  is  simply  to  modify  the  frequency  domain 
equations  derived  in  the  last  chapter  so  that  they  are  written  in  terms  of  the  DFT. 
Since  the  model  is  essentially  a  linear  two  dimensional  convolution  a  certain  amount 
of  zero  padding  is  necessary  to  break  tiie  circularity  of  the  DFT. 

The  method  used  to  find  the  DFT  is  identical  to  that  used  in  section  3.3. 
Acccndingly,  in  this  section,  we  will  tmly  present  the  equations  that  are  different.  In 
order  to  obtain  a  sensible  result,  the  DFT  of  eta,  epsilon,  x(t),  and  the  kernel 
sequence  must  be  evaluated  at  the  same  points  in  the  frequency  domain.  Since  the 
random  variates,  die  time  series,  and  the  kernel  sequence  are  all  of  different  length, 
they  must  each  be  padded  with  zeroes  to  achieve  a  uniform  length  over  which  to 
compute  the  DFT.  Let  the  kernel  sequence  be  N  x  N.  The  time  series  must  then  be 
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at  least  a€  length  2N  making  the  length  of  epsilon  and  eta  3N-1  points.  Further  zero 
padding  may  be  necessaiy  in  ofder  to  efficiently  use  an  FFT  algorithm.  Letting  the 
longest  length  be  n,  the  DPT  d  the  model  can  be  written  as 

Xik)  =  £(*)  +  n(k)  +  ^  £  £(ik)  m-k , )  G\kM  1 )  .  (4.3) 

*,  =  o 


M 

2 


The  underline  signffies  a  zero  padded  sequence.  This  will  produce  a  time  series  of 
length  n.  Only  the  flrst  2N  values  are  retained;  the  remainder  are  either  zero  or 
unwanted  numbers. 

4.3  A  Demonstration  of  the  Inversion  of  a  Bispectrum  with  the 
UBM 

The  software  to  implement  the  inversion  was  written  in  Fortran  77  on  the 
CYBER  830  at  ARL:UT.  It  allows  one  to  specify  the  desired  bispectrum  and  the 
desired  resolution.  Current  limitations  of  the  CYBER  restrict  the  resolution  of  the 
specified  bispectrum  to  a  complex  array  of  128  x  128.  For  convenience,  let  the 
frequency  interval  between  points  be  1  Hz .  This  choice  is  entirely  arbitrary.  Any 
other  number  could  just  as  easily  be  used.  This  restriction  implies  that  the  maximum 
possible  resolution  in  the  bispectrum  is  obtained  when  tlie  isosceles  triangle  in  the 
principal  domain  has  a  base  of  65  points  and  a  height  of  33.  This  in  turn  will 
produce  64^  nonzero  weights. 

Another  factor  is  the  computation  time.  Each  iteration  of  the  specified 
bispectrum  requires  two  128-point  two-dimensional  Fourier  transforms  plus 
additional  computations  on  a  significant  portion  of  the  array.  In  the  light  of  these 
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considerations,  SO  iterations  were  perfnmed  on  the  specified  bispectrum  in  the 
results  presented  below.  This  resulted  in  an  error  of  less  than  1  percent  at  the 
specified  pdnts  of  die  bispectrum. 

4.3.1  The  Results 

For  this  exaniple  the  initial  bispecmim  was  specified  to  be: 

I.  /2)  =  Wi  -  8)  S(f2-  4)  (4.4) 

The  maximunri  resolution  was  used.  A  plot  of  the  resulting  modified  bispectrum  is 
shown  in  fig.  4.3.  This  is  the  bispectrum  that  has  a  bicovariance  which  is  zero  for 
lags  greater  than  63.  The  kernel  sequence  was  found  firom  this  bispectrum.  Using 
the  relation  in  (4.3),  60  128>point  ensembles  were  computed.  The  first  eight 
records  of  the  resulting  time  series  are  shown  in  fig.  4.4.  The  estimated  power 
spectrum  and  the  actual  power  spectrum  are  shown  in  fig.  4.5  and  fig.  4.6.  The 
actual  power  spectrum  was  computed  fiom  the  relation  given  in  section  3.2.3.  As  the 
overlay  plot  in  fig.  4.7  shows,  there  is  very  good  agreement  between  the  estimated 
and  the  actual  power  spectrum.  Since  the  power  spectrum  is  much  easier  to 
compute,  it  was  used  as  an  indicaten’  of  how  many  time  series  were  necessaiy.  Poor 
agreement  here  would  be  indicative  of  even  worse  agreement  in  the  bispectrum.  Fig. 
4.8  shows  an  estimate  of  the  bispectrum  obtained  using  the  Hinich  estimator  [Ifinich, 
1982].  While  it  is  not  helpful  to  overlay  the  initial  bispectram  and  this  estimate,  one 
can  observe  that  they  are  in  good  general  agreement.  By  raising  the  threshold  of  the 
plot  one  can  verify  the  location  of  the  big  peak.  Such  a  plot  is  shown  in  fig.  4.9. 
From  this  plot  it  is  easily  seen  that  the  peak  is  in  the  pre^r  location. 


“Lstimaud  power  spectrum 

Toxtfer  spectrum  computed 
from  the  weights 
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identical  bispectra.  This  bispectrum  is  shown  in  fig.  4.12.  Because  it  is  not 
possible  to  arbitrarily  q)ecify  both  the  bispectrum  and  the  power  spectrum  at  the  same 
time,  it  is  not  possible  to  significantly  change  the  power  spectrum.  One  cannot 
remove  power  at  any  of  the  three  frequencies  associated  with  a  point  in  the 
bispectmm.  Thus  in  essence  all  that  has  been  done  in  the  above  example  is  to  "bury" 
the  peaks  seen  in  the  spectrum  of  fig.  4.11  by  adding  a  very  large  peak  to  the 
spectrum  of  fig.  4.10. 


CHAPTER  5 


SUMMARY 

This  thesis  has  presented  a  method  of  producing  or  estimating  a  stochastic 
time  series  dtat  has  any  desired  bispectral  characteristic.  This  involved  the  use  of  an 
infinite  order  moving  average  noodel  whose  bicovariance  was  obtainable  as  a  finite 
linear  relation  of  the  kernels.  This  expression  was  then  inverted  to  produce  an 
expression  giving  the  kernels  in  terms  of  the  bicovariance.  However,  this  inversion 
is  not  unique.  In  otder  to  identify  a  unique  kernel  function,  an  additional  constraint 
had  to  be  invoked.  This  constraint  was  inqrosed  in  the  form  of  a  relation  between  the 
kernel  function  g(r,$)  and  gfsj*).  The  particular  relation  used  was  g(r,s)  =  g(s,r). 
This  was  an  entirely  arbitrary  choice  used  only  because  of  convenience.  Any  other 
choice  could  have  been  used.  As  was  demonstrated,  a  different  choice  will  not 
change  the  bicovariance,  but  it  will  in  general  change  the  other  polyspectra.  The 
effect  on  the  power  spectram  of  a  different  kernel  relation  is  shown  at  the  conclusion 
of  chapter  4. 

Once  the  kernels  have  been  found  the  problem  is  theoretically  solved. 
Knowing  the  kernels  allows  one  to  compute  the  time  series.  However,  when  one 
attempts  to  implement  this  solution  some  other  problems  arise.  These  problems  are 
the  direct  consequence  of  approximating  infinite  sequences  with  fmite  ones.  In 
general  the  inverse  transform  of  the  bispectrum  is  an  infmite  sequence.  Thus  blindly 
inverse  transforming  the  bispectrum  will  in  general  give  one  an  aliased  version  of  the 
bicovariance.  This  problem  was  circumvented  by  finding  a  bispectrum  which  has  a 
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finite  bicovariance.  This  simply  amounted  to  giving  up  control  over  the  value  of 
some  bispectral  points  in  return  fm  control  over  an  equal  number  of  points  in  the 
bicovariance.  This  was  done  in  such  a  way  so  as  to  have  minimal  effect  on  the 
modulus  of  the  desired  bispectrum.  It  should  be  emphasized  that  if  the  bicovariance 
is  finite  it  can  be  obtained  fiom  the  invose  transform  of  the  bispectrum,  and  in  an 
analogous  manner  the  unaliased  autocoirelation  function  can  be  found  firmn  the  power 
spectrum  provided  the  number  of  points  used  to  represent  the  power  spectrum  is 
equal  to  or  greater  than  the  length  of  the  autocorrelation,  i.e.,  F^l/T  where  F  is  the 
resolution  of  the  power  spectrum  and  T  is  the  length  of  the  autocorrelation.  This 
may  come  as  a  surprise  to  one  thinking  in  terms  of  estimating  spectra  from  a  time 
series.  This  is  the  result  of  confusing  an  estimate  of  the  bicovariance 
(autocorrelation)  with  the  true  bicovariairce  (autocorrelation). 

Having  found  the  kernels  it  turned  out  to  be  easier  to  compute  the  transform 
of  the  weights  and  then  compute  the  DFT  of  the  time  series.  This  is  because  die  DPT 
involves  a  single  integral  (summation)  over  a  product  of  the  transform  of  the 
weights  and  the  Gaussian  variates  as  opposed  to  the  double  integral  (summation)  that 
is  required  in  the  time  domain.  Results  of  a  sample  inversion  were  shown  in  the  case 
where  the  bi^iectrum  consisted  of  a  single  peak.  Other  inversions  were  also  done 
yielding  equally  good  results.  Generally  speaking,  increasing  the  variance  increases 
the  number  of  time  series  needed  to  reproduce  the  bispectrum.  This  is  simply  due  to 
the  fact  that  ^silon  and  eta  are  produced  by  pseudo-randmn  number  genmtors.  The 
resulting  time  series  has  the  propo’  bispectrum  to  the  extent  that  the  random  variates 
approach  the  Gaussian,  zero  mean  assumption.  A  higher  variance  simply  means  that 
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nme  nmdom  variate  records  are  needed  befcne  the  estimated  third,  fifth,  and  higher 
order  mooients  are  actually  zero. 
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APPENDIX  A 

‘  COMPUTING  THE  BICOVARUNCE  WITH  THE  TWO- 

DIMENSIONAL  FFT 

■  Since  the  bicovariance  is  derined  over  the  interval  (-N,  N)  and  the 

'  bispectrum  is  defined  over  frequencies  [-1/2,1/21,  they  must  both  be 

appropriately  shifted  before  being  transformed  with  the  DPT  so  that  B(0,0)  or 

( 

C(0,0)  sum  without  being  phase  shifted  and  not  B(- 1/2,- 1/2)  or  C(-N,-N).  The 
mechanics  of  this  translation  are  shown  in  fig.  A-1.  By  translating  the  bispectrum 
in  this  way,  the  negative  fiequency  term  {exp(-j2ic(kili-i-k2l2)/N))  of  the  DPT  is 
being  replaced  by  the  equivalent  term  {exp(j2jc((N-ki)li+(N-k2)l2)/N)).  The 
inverse  transform  of  this  array  is  an  analogously  shifted  version  of  the 
bicovariance.  The  mechanics  of  the  translation  of  the  bicovariance  is  nearly 
identical.  The  only  difference  is  that  the  bicovariance  must  be  padded 
two-dimensionally  with  a  zero  in  cnder  to  make  its  dimensions  even  so  that  the 
FFT  may  be  used. 
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APPENDIX  B 


TIME  SERIES  GENERATION  SOFTWARE 


C  WRITTEN  BY  PETER  ALLISON. 

C  JANUARY  1988 

C - 

C  THISmOGRAMCXJMPinES  A  TIME  SERIES  WITH  A  GIVEN  BISPECTRAL 

C  CHARACTERISTIC  BY 

C  1)  COMPUTING  THE  WEIGHTS  IN  THE  FREQ  DOMAIN  FROM  THE  GIVEN 

C  BISPECTRUM 

C  2)  COMPUTING  THE  FT  OT  THE  DESIRED  SIGNAL  IN  THE  FD  (FROM  THE 

C  WeGHTS) 

C  3)  INVERSE  DF  TRANSFORMING  THE  RESULT  TO  PRODUCE  THE  DESIRED 
C  SIGNAL 

C 


PROGRAM  BTSFTN4 

INTEGER  IRCJfllECJWK(9l8).TESTJ42>UC2-i.TS  JfV 
INTEGER  NRCERR 

REAL  XT(128)XS(1024)XP(1024)ETA(200).GNU(200) 

REAL  SIG^CALE>!AXAlAXS>flNSJ>SE(128)4^T(512)JlWK(918) 
REAL  TIC 

DOUBLE  PREQSION  DSEED 

COMPLEX  CX(128).CE(128),CN(128).SUM.C(128.128).CG(I28.128) 
COMPLEX  CXTS(128),CWK(128) 

OPEN(UN1T«94TLE=’WEIGHTO 

DATA  PSE/128*07 
DATA  MINS>IAXS>IAX/3*07 
DATA  C.CG/32768*(0.,0.y 
DATA  CXTS/I28*(0..0.)/ 

DATA  ERR/1/ 

REWIND(9) 


C  THERE  IS  AN  OPTION  TO  EITHER  READ  THE  BICOVARIANCE  IN  FROM 

C  A  FILE  (WEIGHTQ  OR  TO  READ  THE  BISPECTRUM  IN  FROM 

C  FILE  DATA. 


4  PRINT*,'INPUT1  TO  GENERATE  THE  MODEL  WEIGHTS,  INPUT  2' 
PRINT*;  TO  READ  THEM  IN  FROM  FILE  WEIGHTC 
READ*,TEST 
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WUNT*.’  • 

LTS«128 

IFCrEST£Q.l)THEN 

2  PIUNT*,TNPUT  THE  NUMBER  OF  WEIGHTS  TO  USE* 

PRINT*,*  (MAX  64)* 

READ*J< 

IF(N.GT.64)*IHEN 

PRINT*,'  NUMBER  IS  GREATER  THAN  64* 

GOT02 
ENDIP 
PRINT*,*  * 

N2«2*N 

XCALLREWD 

ENDIF 

PRINT*,*INPUT  THE  NUMBER  OF  RECORDS  TO  GENERATE* 

READ*,NREC 

PRINT*,*  * 

HUNT^.TNPUT  SIGMA* 

READ*,SIG 
PRINT*,*  * 

IF(TESTiQ  Jl)  THEN 
READ(9,*)N 
N2»2*N 

PRINT*,*ENTER  0  TO  DISCARD  EXISTING  T.S.* 
READ*JNRC 

IF(NRCJEQ.O)  CALL  REWD 
CALL  GENERG(C.CG,N23IG  JIWKJWK,CWK,TEST) 
ELSE 

1  CALL  BISPIN(C,CG  J«,SIGJWK41WK,CWK,TEST,ERR) 

IF(ERR£Q.l)THEN 

PRIN'P.'ENTER  A  BETTER  VALUE  FOR  N* 

READ*J4 
PRINT*,*  * 

N2*2*N 

GOTOl 

ENDIF 

ENDIF 

CLOSE(9) 

OPEN(UNIT»44TLEa'BTSOUT') 

OPEN(UNrr-54=ILE«TIMSER*) 

OPEN(UNrr-6,FILE*TIMSERI’) 

REWIND(4) 

WRrrE(4,*)  ’  NUMBER  OF  WEIGHTS  USED  IS  'J<i 
WRrTE(4,*)  *  NUMBER  OF  RECORDS  COMPUTED  IS  *J^REC 
WRITE(4,*)  *  VARIANCE  OF  ETA  AND  EPSILON  IS  '.SIG**2 

C 
C 
C 
C 


BISPECTRUM  HAS  BEEN  LOADED  INTO  ARRAY  C. 

THE  UNAUASED  BICOVARIANCE  HAS  BEEN  COMPUTED. 
THE  MODIFIED  BISPECTRUM  IS  NOW  IN  ARRAY  C. 


C  FOURIER  IKANSFORMCX’ WEIGHTS  IS  IN  ARRAY  CG 
C  NOW  GENERATE  THE  TIME  SERIES  IN  THE  FREQUENCY  DOMAIN 
C - 


PRINT*.'GENERATE  RANDOM  VARIATES' 

C»^(UNIT»3JTLE-'DSEED') 

REWIND(UNrr«3) 

READ(33)  DSEED 
3  FCXIMAT(D17.10) 

IF(N.LT.42)  THEN 
FPAD-0 
NV»3*N-1 
ELSE 

FPADsl 

NV=LTS-1 

ENDIF 

C  LOOP  110  IS  DONE  ONCE  FOR  EACH  RECORD  GENERATED 

DO  110IRC=1J^C 

CALL  GGNML(DSEEDNV.GNU) 

CALL  GGNML(DSEEDNV^A) 

DO120I=l>rV 

gnu(i><;nu(i)*sig 

ETA(I)i^A(I)*SIG 

CE(I)=CMPLX(ETA(D) 

CN(D=CMPLX(GNU(D) 

120  CONTINUE 

DO130I=NV+l^TS 

CE(I>nCMHJC(0.) 

CN(I)=CMHLX(0.) 

130  CCWTTNUE 

CALL  FFT2C(CE.7  JWK) 

CALL  FFT2C(CN,7  JWK) 

C  THE  DFT  OF  THE  PADDED  VARIATES  IS  COMPUTED  IN  120  AND  130 

DO140Ka0a-TS-l 

CX(K+1)=CMPLX(0) 

SUM=CMPLX(0) 

DO150L=0^TS-l 

M=K-L 

IF(MiT.0)  M=LTS+M 

SUM=SUM4CE(L+1)*CN(M+1)*C0NJG(CG(L+1.M+1)) 

150  CC»mNUE 

CX(K+1)=SUM/LTS  +  CE(K+1)  +  CN(K+1) 

CXTS(K+ 1)=CXTS(K+1)+CX(K+1) 

CX(K+1)=C(»GG(CX(K+1)) 

CXTS(K+l)=CONJG(CXTS(K+l)) 

PSE(K+l>=REAL(CX(K+l)*CONJG(CX(K+l)))A.TS+PSE(K+l) 

140  CCmiNUE 

13  PORMATCSUMC,I2.>'.4(E16.7)) 
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C  TAKE  THE  INVERSE  TRANSFORM  OF  THE  TIME  SERIES 

CALL  FFT2C(CX.7  JWK) 

WRrrE(4  .•)  ••••  XT(D  ••••  cx(D  cxa+N)  •••••■ 

DO160I-1J-TS 

XT(I)a«EAL(CX(D)/LTS 
IF(IRCLT.9)  XS((IRC-1)*LTS+I)»XT(D 
IF(MAXSI.T.XT(I))  MAXS»XT(I) 

IF(MINS.GT.XT(I))  MINS»XT(I) 

WRnE(4.11)  LXT(I).CX(I)^ 

IF(FPAD£Q.O)THEN 
WRnE(6.9)XT(D 
IF(LGT^*N)  XT(I)=0. 

ELSE 

WRrTE(6^)XT(D 

IF(LGTLTS-N)XT(I)=0. 

ENDIF 

WRnE(5^)XT(D 
160  CONTINUE 

HUNT^.TINISHED  THE  ‘JRC  .TH  TIME  SERIES' 

1 10  CONTINUE 

1 1  FORMATC  XTC43.VJ8.3.*  CX(I)=’;2(E12.4).'  K  =  ’43) 
CLOSE(5) 

CLOSE(6) 

REWIND(UNIT*3) 

WRrrE(3.3)  DSEED 
CLOSE(3) 

C  COMPUTE  THE  BISPECTRUM 

CALL  BISPGEN(CG.C.CXTS4*TS,SIGJ«EC) 

CALL  FFnC(CXTS.7  JWK) 

PRnsTP.TINISHED  COMPUTING  THE  BISPECTRUM’ 

OPEN(UNrr=7FlLE*’TIMSER2') 

DO  170  I*14-TS 

CXTS(I)*CXTS(I)/(LTS'»NREC) 

WRnE(7^)  REAL(CXTS(I)) 

170  CONTINUE 
CLOSE(7) 


PRINT*,’  BEGIN  POWER  SPECTRUM  COMPUTATION’ 
WRITE(4.*)  ■***••••••  POWER  SreCTRUM  •••••••••■ 

DO  200  I=lJLTS/2+l 
PSE(I)=PSE(D/NREC 
WRITE(4.8)  IJ«E(I) 

IF(PSE(I).GT.MAX)  MAX=PSE(I) 

200  CONTINUE 

DO  210 1=1,1024 


XP(I>=I 

210  CONTINUE 


c  MjOt  the  time  series  and  the  bispectrum 

CALL  PLTLFN(L"PLOTrj-"1234567890M0000) 

CALL  PLTORG(I..I.) 

CALL  PLTAXIS(0..0.,5..0..0..65.^.X"FREQUENCY".-  J  04) 
CALLPLTAXIS(0..0.4..90..0..1.,.2^"EST  PS  MO.l) 
CALLPLTDATA(XPJ>SE.65.0,0.0.,65y5..0..MAX/5...08) 

CALL  PLTEND(8.5.1I.O) 

CALL  P0WSPEC(CG.SIG4>ST,MAX) 

CALL  PLTLFN(LTLOTI"J-"I234567890", 10000) 

CALL  PLTORG(l..l.) 

CALL  PLTAXIS(0..0.,5..0..0.,65..2.4-"FREQUENCY", -104) 

CALL  PLTAXIS(0.,0.4..90..0..1...2^"COMP  PS  MO.l) 
CALLPLTDATA(XPJ>ST.65,0,0,0..65y5..0..MAX/5.,.08) 
CALLPLTEND(84,11.0) 

Dr  (NREC.GT.7)  THEN 

PRINT*, MAX  IS  'MAXS;  MIN  IS  'MINS 
PRINT*, INPUT  STARTING  VALUE  FOR  PLOT 
READ*M1NS 

PRINT* .’INPUT  THE  TIC  INTERVAL’ 

READ*,TIC 

IF(MrNS.GT.0)MINS=O. 

SCALE=(MAXS-MINS)/5 

CALL  PLTLFN(L  ”PLOTrjL’'1234567890M0000) 

CALL  PLTORG(l..l.) 

CALL  PLTAXIS(0..0.4..0.,0.,1024.,20.,L"’nME",- 10,5) 
CALLPLTAXIS(0.,0.4..90..MINS,MAXS,'nca."AMPLITUDEMO,2) 
CALLPLTDATA(XP,XS,1024,0.0,0.,1024/5.,MINS,SCALE,.08) 
CALL  PLTEND(84.11.0) 

ENDIF 

8  FORMAT(I3.2(E17.8)) 

9  FORMAT(E17.8) 

OPEN(UNIT=  18  JTLEs'BISPDAT) 

REWIND(18) 

READ(18,*)  NRC 

NREC=NRC+NREC 

REWIND(18) 

WRITE(18,*)  NREC 
CLOSE(18) 

CALL  BSDATFL(LTS4JREC) 

END 


nnoonn 
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C  •••*•••  SUBROUTINES 


SUBROUTINE  BISP1N(C.CGJ4.SIG4WKJIWK.CWK.TEST.ERR) 


THE  BISPECTRUM  IS  INPUT  FROM  FILE  DATA  THEN  THE  UNALIASED 
BICOVARIANCE  IS  FOUND  BY  ITERATION.  THE  BISPECTRUM. 
THE  TRANSFORM  OF  THE  WEIGHTS.  AND  THE  BICOVARIANCE 
ALL  ARE  WRITTEN  OUT  TO  FILES. 


REAL  CST(64).SIGJIWK(918) 

INTEGER  L1(64)J-2(64)4=J42.ERR.TEST.IWK(918) 

COMPLEX  C(128.128).CG(128.128)3.CWK(128) 

OreNCUNTTsi  1  ^TL^ETOUT) 

OPEN(UNIT=2JTLE=‘DATA') 

OPEN(UNIT*15JFILE=PBISP0 

REWIND(UNIT=2) 

REWIND(UNIT=1) 

REWIND(15) 

WRITEdS.*)  N/2+1 
WRITEdS.*)  N/4+1 
N2=N/2 

PRINT*  .INPUT  NUMBER  OF  NONZERO  BISPECTRUM  POINTS’ 

READ(2.4)  NP 

4  F0RMAT(I3) 

PRINT».'NUMBER  OF  NONZERO  BISP.  POINTS  IS  ’J4P 
F=N+1 

HUNT*  .INPUT  BISPECTRUM  FROM  FILE  DATA' 

DO  5  1=1 

READ(2.7)L1(I).L2(D.CST{I) 

PRINT*  .’BCJL1(I)JL2(I).’)=  '.CST(D 

IFaia)GTJ9/2.0Ri2(I).GTJi/4.0Ril(I)+L2(I).GTJ9/2)THEN 
PRINT*.’POINT  OUTSIDE  THE  PRIN.  DOMAIH 
ERRsI 
GOTO  99 
ENDIF 

5  CONTINUE 

7  FORMAT(I3.I2E4.1) 

PRINT*.'INPUT  #  OF  ITERATIONS' 

READ*JC 

WRITEd.*)  •••••♦••**•*•  K  *  -jc.’**************' 

C  EACH  LOOP  30  REPRESENTS  ON  ITERA’nON  OF  THE  BISPECTRUM 


DO  30  M=1  JC 
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CALLNEXISTB(CJ^ 

DO  15 1=1  JiP 
B=CMPLX(CST(I)) 

C(L1(I)+1J.2(I)+1)=B 

C(L2(I)fli-l(I)+l)=B 

C(Ll(I)+L2(I>fl.-Ll(I>fF)=B 

IF(L2(I)J1E.0)  C(L1(I)+L2(I)+1.'L2(I)+F)=B 

IF(L2(DJ«.0)  C(-L2(D+F.-L1(D+F)=B 

IF(L2(I)J^.O)  C(-L1(I)+F.-L2(I)+F)=B 

C(-Ll(I>L2(I>fF^2(D+l)=B 

C(L2(I>f  1JF-L1(1>L2(I))=B 

C(Ll(Dtlf-Ll(I)-L2(I)>*B 

C(F-Ll(I)-L2(I)JLia)+l>»B 

IF(L2(DmO)  C(F-L2(I).L2(I>t-Ll(D+l)=B 

C(F-Ll(I)a-l(I>+L2(I)+l)=B 

15  CONTINUE 

CALL  FFT3D{C.128.128J4J<.1.-1.IWKJ1WK.CWK) 

IF(M.EQ.K)  THEN 
DO  10I=1J42 
DO  10J=1J42 

10  WRrrE(i.i)i-U-i.caJ)>i 

ENDIF 

C  BEFORE  TRANSFORMING  SET  TO  0  VALUES  OF  BICOV.  THAT  DONT  EXIST. 
IF(MJ4E.K)  CALL  NEXISTC(C J4) 

C  IF  ms  THE  LAST  LOOP  THEN  USE  THIS  BlCOVARIANCE  TO  COMPUTE  G 

IF(M.EQ.K)  CALL  GENERG(C.CGJ^,SIG J1WK4WK.CWK.TEST) 

IF(M.EQ.K)  THEN 
DO  50 1=1.128 
DO  501=1,128 

WRITE(1.3)  I-U-l.CGaJ) 

50  CONTINUE 

ENDIF 

CALL  FFT3D(C.128,128J<J^,l,l.IWKll'WK,CWK) 

IF(M.EQ.K)  THEN 
00201=11^2/2+1 
DO20I=lJ42+l 

IF(1.LE.I.ANDJ+IXE.N2)  THEN 

WRITE(15.*)  SQRT(REAL(CaJ)*CONlG(Cai)))) 

ELSE 

WRrrE(i5.*)0. 

ENDIF 

20  CONTINUE 

WRrrE(l,2)  I-U-I.ca4)>l 
ENDIF 

1  FORMATC  C(M2;.M2.')  =  '.2(E16.6),'  K=  ’.12) 

2  FORMATC  BC.I2.’.’12,')  =  '.2{E16.6),'  K=  ’,12) 

3  FORMATC  GC43.’.’13.’)  » ’.2(E16.6)) 

PRINT*.’F1NISHED  THE  ’M.TH  ITERA-nON  OF  THE  BISPECTRUM’ 

30  CONTINUE 
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cpB^ 

99  CONTINUE 
CLOSEd) 
CL0SE(2) 
CLOSEdS) 
RETURN 
END 


SUBROUTINE  GENERG(C.CG  J<,SIG  JIWKJWK.CWK.TEST) 


C - 

C  TfflS  SUB  COMPUTES  THE  KERNELS  FROM  THE  BICOVARIANCE 
C  IT  USES  A  SYMMETRICAL  RELATICW  BETWEEN  TOE  WEIGHTS. 

C  HOWEVER  ANY  RELATION  COULD  BE  USED. 

C - 

COMPLEX  CG(128.I28).C(128.128).CWK(128) 

REALRWK(918).SIG 

INTEGER  NJWK(918).TEST 

OPEN(UNIT=8J!ILE=WEIGHTS’) 

REWIND(UN1T=8) 

C  EITHER  THE  WEIGHTS  ARE  READ  IN  OR  THE  BICOVARIANCE 

C  COMPUTED  IN  SUB.  BISPIN  IS  WRITTEN  OUT  TO  A  FILE 

IFfTESTXQ.l)  THEN 
WRITE(8.*)  N/2 
WRITE(8.*)  N/2 
DO30I=lJM/2 
DO30J=lJ^/2 

WRITE(8.*)  REAUCa  J)) 

30  CONTINUE 

ENDIF 

IF(TEST.EQ.2)  THEN 
PRINT*.TNPUTINGG’ 

D0  5I=lJI/2 
D0  5J*lJ</2 

READ(9.*)  VALUE 
CaJ)=CMPLX(VALUE) 

5  CC8STINUE 

ENDIF 

C  FIND  THE  WEIGHTS 

DO  10  I=2JM/2 

CG(1J)=.25*C(14)*(17S1G**4) 

CGa.l)=CG(lJ) 

10  CONTINUE 

CG(l.l)=(W6.)*C(l.l)*(iySIG**4) 

DO  20 1»2  J</2 
DO20J=2J^/2 

CGaj>=.5*C(IJ)*(l.^IG**4) 

20  CONTINUE 

PRINT*.'COMPUTING  THE  TRANSFORM  OF  THE  WEIGHTS’ 

CALL  FFT3D(CG.128.128.128.128.1.1.IWKJIWK,CWK) 

3  FORMATC  GC.I3.’.’J3.7  =  '.2(E16.6)) 

DO  501*1.128 
DO50J»1.128 

50  WRrTE(1.3)  I-l  J-l.CGOJ) 
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CL0SE(8) 

RETURN 

END 


Oonoo 
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SUBROUTINE  BISPGEN(CX}.C.CXTSJ.TS^IGJ4R£C) 


THIS  SUB  COMPUTES  THE  BISPECIRUM  FROM  THE  WEIGHTS  AND  IT 
COMPUTES  TTE  BISPECTRUM  OP  THE  AVERAGED  TIME  SERIES. 
THIS  IS  A  HIGH  VARIANCE  ESTIMATE  OF  THE  BISreCTRUM. 


COMPLEX  CXTS(128).CG(128.128)4»MET.C(128.128) 

INTEGER  LTS  Xl  JC>IJJ  JfllEC 

REAL  TEMP^IG 

PRINT*.'COMPUTING  BISPECTRUM* 

OPEN(UNIT*16JTLE.PBISPr) 

OPEN(UNIT=I7JILE*PBISPE') 

REWIND(16) 

REWIND(17) 

CALL2ER02(C.128) 

K.LTS/4 

DO  10 1-OXTS/2 
IFa.LE.K)THEN 
M-I 
ELSE 

M=iLTS;/^-I 

ENDEF 

DO20J=0M 

IF(JJ^.0)THEN 

L2»128-J-t-l 

ELSE 

L2*l 

ENDIF 

IF(I  J^.O)  THEN 
Ll»128-1+1 
ELSE 
Llal 
ENDIF 
U=I+J+1 

PMETCGa+1  J+l)+CONJG(CG(Ll  JJ)>+CONJG(CG(U  1)) 

PMET*PMET+CX>NJG(CG(L24;))+C0NJG(CG(UL.2))4CG(J+1.I+1) 

WRrrE(16*)PMET 

PMET«CXTS(I+1)*CXTS(J+1)*C0NJG(CXTS(IJ)) 
WRITE(17*)PMET 
20  CONTINUE 

DO30J»M-t-l.K 

PMET»CMPLX(0.) 

WRITE(17*)PMET 
WRrrE(16*)PMET 
30  C(»mNUE 

10  CONTINUE 


no 
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THE  REMAINING  COIS  SIMPLY  WRITES  THE  DATA  OUT  IN  A  FORM  THAT 
CAN  BE  PLOTTED. 

REWIND(16) 

REWIND(17) 

DO40I-0J-TS/2 
DO40J-0JC 

READ(16*)Ca+IJ+l) 

40  CX»mNUE 
REWIND(16) 

WRrre(16*)K+l 
WRITE(16  .*)  LTS/2+1 
DO  50  Is«0JLTS/2 
DOSOJ«KA-1 
PMET-Ca+U+I) 

TEMPKSIG*M)*S<»TCREAL(n>IET*a»IJG(W4ET))) 

WRrrE(I6*)'I1EMP 
50  OWTINUE 
DO60I*0JLTS/2 
DO60J«0JC 

READ(17*)Ca+U+l) 

60  CWTINUE 
REW1ND(17) 

WRITEd?,*)  K+1 
WRTIEd?  ,•)  LTS/2+1 
DO  70 1«0J-TS/2 
DO70J«K,0.-l 
PMET«C(I+U+l) 

TEMP»REAL(PMET*C(»OG(PMEr)yNREC**3 
WRITE(17,*)TEMP 
70  CONTINUE 
CLOSE(16) 

CLOSE(17) 

RETURN 

END 


on 
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SUBROUTINE  REWD 
C»>EN(UNrr«18JIL£='BISPDAT) 
OPEN(UNIT«7JTLE»TIMSER2') 
REWIND(5) 

REWIND(6) 

REWIND(7) 

R£WIND(18) 

WRITE(18  .•)  0 
CL0SE(18) 

CL0SE(7) 

RETURN 

END 


SUBROUTINE  POWSPEC(CG.SIGJPST>lAX) 

COMPLEX  CG(I28.I28) 

REAL  PST(128).SIG;SUM>IAX 
MAX=0 
DO  10 1=0,127 
SUM=0 
DO  201=0,127 

M-I-J 

IF(MiT.O)  M=128+M 

SUM=SUM+REAL(CG(J+1>1+1)*C0NJG(CG(J+1>1+1))) 
20  CX»mNUE 

PSTa+l)=2*SIG**2+SUM*(SIG**4) 

IF(PSTa+l).GT.MAX)  MAX=PSTa+l) 

10  CONTINUE 
RETURN 
END 


SUBROUTINE  ZERO(A  J4D) 

COMPLEX  A(128) 

DO  10 1=1 

A(D=CMnJC(0) 

10  CONTINUE 
RETURN 
END 

SUBROUTINE  BSDATFL(N  JIREC) 

TraS  SUB  CREATES  THE  NECESSARY  FILES  WHICH  CAN  BE  READ  BY  THE 
BISPECIRUM  ESTIMATION  ALGORITHM 

INTEGER  NUMJ^  JTOEC 
NUM=N*NREC 

OPEN(UNIT=10JFILE=REFINF) 

OPEN(UNIT=  1 1  Jra-EsEFTINF) 

OPEN(UNIT=12JTLE='PWRGINF7 
OPEN(UNrr^l3JTLE=TWRAINF^ 

OreN(UNIT=  I44TLE=’BISPINF^ 


REWIND(10) 

REWINDCH) 

REWIND(12) 

REWIND(13) 

REWIND(14) 

WRITE(10.*)  "TIMSER” 
WRITE(10 .*)  "TSPIO* 
WRITEdO .•) 

WRITEdO.*)  ‘l.l.l’ 

WRITE(10.*)  N 
WRITEdO  •)  "REAL" 
WRITEdO,*)  "REAL" 
WRnE(lO,*)  T' 

WRITEdO,*)  "(E17^)" 
WRTIEdO,*)  TT 
WRITEdl,*)  "TSPIO" 
WRnE(ll,*)  *T=FT~ 
WRnE(ii,*)  T;j<REC,',r 

WRITEdl,*)  N _ 

WRITE(12,*)  "FFT" 
WRnE(12,*)  "POWSPC" 
WRrrE(l2,*)  T.*J«EC,*,1* 
WRnE(l2,*)  2*N 
WRITEd3,*)  "POWSPC" 
WRnEd3,*)'"NEWPWR" 
WRnEd3,*)  T,'JIREC,',r 
WRrrEd3,*)  N/2+1 
WRrrEd3,*)  NREC 
CLOSEd3) 

OPEN(UNIT»144=ILE=BISPINF) 
WRITE(14,*)"’SPIO" 
WRnE(l4.*)"FFT" 
WRnE(l4,*)  T,'J«EC,’,r 
WRnEd4.*)  T,',2*N,',r 
WRITE(14.*)"Y" 
WRITE(14,*)"SPIO'" 
WRITE(14,*)"NEWPWR" 
WRITE(14.*)  T,l,r 
WRrrE(i4,*)  T;j4/2+i;,r 
WRnE(14,*)  N 
WRITE(14.*)  •2' 
WRnE(l4,*)"N" 

WRnE{l4,*)  NREC 

WRITE(14,*)"BISPRr’ 

WRITE(14,*)"PWAVG" 

WRnE(14,*)"BIRAW" 

WRnE(14,*)™CHI" 

WRITE(14,*)"N’- 

WRITE(14,*)  "BISUM*" 

CLOSE(10) 

CLOSE(ll) 

CLOSE(12) 

CLOSE(14) 

RETURN 
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END 


10 


20 


SUBROUTINE  NEXISTC(CJ^ 
COMK^  C(128.128);£ 
INTEGER  N 
Zf=CMPLX(Q) 

DO  5 1=1 
DO  5  J=1JN 

CaJ)=CMPLX(REAL(C(U))) 
CONTINUE 
DO  10 1=2^/2 
DO  10  J«TO+1 
C(N/2+U>=Z 
C(JJl/2+I)=Z 
CONTINUE 
DO20I=lJ< 

C(N/2+lJ)=Z 

caj^/2+i>z 

CONTINUE 

RETURN 

END 


SUBROUTINE  NEX1STB(CJ^ 

INTEGER  NCT 
COMPLEX  C(128.128) 

DO  10 1=1 
DO  10J=1J^^ 

CaJ)=CMPLX(0) 

10  CONTINUE 

C  ZERO  OUT  THE  BISP  THAT  DOESNT  EXIST-  IF  NYQUIST  IS  TO  BE  SATISFIED 


NCT=2 

DO20J«N/2-lJ^Al+l.-l 

DO30I=NCTJ 

Ca+U+1)=CMPLX(0) 
C(J+U+1)=CMHA(0) 
30  CONTINUE 

NCT=NCT+1 
20  CONTINUE 
NCT=N-2 

DO40J=N/2+lJM-2 
DO50I=N/2+l^CT 
C(I+U+1)=CMPLX(0) 
C(J+U+1)=CMPLX(0) 
50  CCKTINUE 

NCT=NCT-1 
40  CONTINUE 
RETURN 
END 
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SUBROUTINE  ZERO2(B>0 
INTEGER M 
COMPLEX  B(128.128) 

DO  10I-1J4 
DOIOJ-IX 

B(IJ)=CMPLX(0) 

10  CONTINUE 
RETURN 
END 
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